Python tutorial

On this Page we give a brief introduction to python, for people who want to be users of python ESPResSo or ESPResSo++. The goal of the tutorial is writing a short program that measures the escape time distribution of a Langevin particle from a potential well.

= Introduction = In this tutorial the syntax, basic control structures, datatypes and a bit of plotting with matplotlib and numerics with numpy can be learned.

Python's is a (mostly) procedural/object oriented language and is easy to learn for anybody who is familiar with languages like C, C++ and Fortran. Lots of material is found online and to a large extend of learning about the powerful features of Python libraries is actually searching the web for how to use them. This is why we do not give detailed explanations, but hope that the reader learns easily by example.

= The Langevin Equation = A bit of motivation towards that.

= A route towards the Langevin particle = Here we suggest a sequence of steps to obtain a program that integrates the Langevin equation of motion and measures the mean escape time from a potential well.

Define potential and force
First you need to define a harmonic potential such as $$V(x) = V_0x^2$$ and a corresponding force $$F(x) = -2V_0x$$. You can start with an empty file and implement a function like:

 def V(x): return V0*x**2 


 * Implement the functions V(x) and F(x)

As a general remark, note that in the following you can set the mass of the particle to unity.

Plot harmonic potential and force
In order to plot the potential you can use for instance the pylab module. We also use the numpy module here to generate an array of positions for the x-axis of the plot.

 import pylab import numpy as np 

A simple plot of the potential and the force can be obtained by calling the function below.

 def plot_VF: x = np.linspace(-5,5,101) pylab.plot(x,V(x)) pylab.plot(x,F(x)) pylab.xlabel('xlabel') pylab.ylabel('ylabel') pylab.title('title') pylab.grid(True) pylab.show 


 * Add and call the plot function in the script

Implement integrator
This will be the inner integration loop of the algorithm which contains the integrator. It is sufficient to use a simple integrator according to the Euler scheme:

$$ x(t+\Delta t) = x(t) + v(t)\Delta t$$

$$ v(t+\Delta t) = v(t) + a(t)\Delta t$$

 def integrate(x,v,nsteps,dt): ...   return (x,v)  that takes the current position x, velocity v and the number of integration steps as well as the the time step as arguments. After integrating for nsteps it should return the current position and velocity of the particle.
 * Write the integration function

Outer Integration Loop
It is comfortable to use an outer integration loop which contains both the integrator and a part that calculates and stores observables. This loop performs the whole simulation, so just call it once at the end of your script.

 def run(nblocks): ... 
 * Create a function
 * Set the initial conditions
 * Prepare the lists for observables that you will use later in the analysis
 * Do the integration
 * Store the time t, position x and velocity v to lists

Use lists to store the observables. Create a list that contains the initial value  pos = [Initial_x_value]  and use the command  pos.append(Next_x_value)  to add more entries. You can now run the whole simulation by calling

run(nblocks)

with the desired number of outer cycles.

Your code shoule now be structured like this:

 def run(nblocks): set initial conditions prepare lists for time, positions and velocities for i in range(nblocks): integrate store t, x and v to lists 

Plot trajectory in time and phase space

 * Plot the trajectory of $$(t,x(t))$$ and $$(t,v(t))$$ in time and in phase space $$(x(t),v(t))$$.

You can either write a function to do this or add the code directly after the outer loop in the run function.

 def run(nblocks): set initial conditions prepare lists for positions and velocities for i in range(nblocks): integrate store t, x and v to lists plot trajectories </tt>

Calculate and plot kinetic and potential energy
A good position to implement the functions is where x and v are stored.
 * Write functions that calculate kinetic and potential energy
 * Store $$E_{kin}$$ and $$E_{pot}$$ to lists.


 * Plot the energies

 def run(nblocks): set initial conditions prepare lists for positions and velocities for i in range(nblocks): integrate store t, x and v to lists calculate and store potential and kinetic energy to list plot trajectories plot energies </tt>

Add a command line parser

 * Add a command line parser for nsteps, nblocks, $$V_0$$ and $$\Delta t$$

You might want to change some settings of your simulation when calling it from the command line instead of modifying your code for each small change. Therefore python comes with a command line parser called optparse.

from optparse import OptionParser

The initialization is done by creating an object "parser".

parser = OptionParser

You can add entries for which the parser will look in the given arguments by using

parser.add_option("-V", "--potential", dest="V0", default=1)

while the parsing itself is done with the following statement:

(options, args) = parser.parse_args

Use

 options.V0 </tt>

to access the variable that was parsed from command line.

Langevin particle
In the following you will modify the program for the harmonic oscillator in order to simulate a particle driven by random forces in a potential with finite height.

Define the stochastic Langevin force
As a first step introduce the Langevin dynamics to the algorithm by adding an additional force term $$F_{Langevin} = -\gamma v + \xi \sqrt{\frac{2\gamma k_BT}{m \Delta t}} $$.

 def F_langevin(v,gamma,T,dt): ... </tt>
 * Add the Langevin force term

Here $$\xi$$ is a gaussian distributed random number. Note that as in the previous part you can set mass $$m$$ and Boltzmann constant $$k_B$$ to unity. Since we are using the force in a discretized simulation we need to divide by the square root of the time step to correct for the collision frequency.

In order to obtain $$\xi$$ there is a module called "random" that is able to generate this type of random number via

 import random xi = random.gauss(mean, width) </tt>

Check temperature (histogram of position)

 * Modify the command line parser to parse also the temperature of the simulation and use this temperature in the Langevin force term.


 * Plot a histogram of the position of the particle

Here you can use the hist( list, bins) function from pylab.  pylab.hist(list, bins) </tt>

Change potential and forces
Here we want to simulate the particle in a potential with finite depth $$V(x) = - V_0 \exp(-\frac{x^2}{\sigma^2})$$.


 * Modify V(x) and F(x) according to the new potential.

Add escape check
Depending on the Temperature that we use for the driving Langevin term the particle might able to escape from the potential dip. This is the case we want to investigate further in the following steps, so choose a temperature that is high enough for the particle to escape. Look at the histogram form the previous section to see what happened.


 * Choose a border that the particle has to pass until it "escaped" from the potential (e.g. twice the width of the Gaussian) and add a query that checks this every integration step.


 * Modify the function "integrate" such that it immediately returns the time needed to escape (and stops to integrate further).

 if particle escaped return escape time </tt>

Since you are no longer interested in the positions and velocities at the end of the simulation remove the lists and the trajectory plots (if not done already...). The outer integration loop is now used to repeat the experiment and calculate the average. Therefore make sure that the setting of the initial conditions do now take place within this loop and also choose dt and nsteps such that the particle does escape.

Measure mean escape time in outer loop

 * Store the escape time of each simulation for the inner loop in a list
 * Calculate the mean and the standard deviation

Here it is very useful, that there are some build-in functions for the type list

 sum(list) # returns the sum of all elements len(list) # returns the length of the list </tt>


 * Print the results at the end of the simulation
 * Repeat the simulation for different temperatures and observe how the escape time scales

Improve code structure

 * Encapsulate the Langevin Sim in a class