A MATLAB program to calculate translational and rotational diffusion coefficients of a single particle

Computer Physics Communications 182 (2011) 400–408

Contents lists available at ScienceDirect

Computer Physics Communications www.elsevier.com/locate/cpc

A MATLAB program to calculate translational and rotational diffusion coefficients of a single particle ✩ Mohammad A. Charsooghi ∗ , Ehsan A. Akhlaghi, Sharareh Tavaddod, H.R. Khalesifard Department of Physics, Institute for Advanced Studies in Basic Sciences (IASBS), Zanjan 45137-66731, Iran

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 20 May 2010 Received in revised form 8 September 2010 Accepted 22 September 2010 Available online 29 September 2010 Keywords: Translational diffusion Rotational diffusion Brownian motion Single particle Diffusion coefficients MATLAB program

We developed a graphical user interface, MATLAB based program to calculate the translational diffusion coefficients in three dimensions for a single diffusing particle, suspended inside a fluid. When the particles are not spherical, in addition to their translational motion also a rotational freedom is considered for them and in addition to the previous translational diffusion coefficients a planar rotational diffusion coefficient can be calculated in this program. Time averaging and ensemble averaging over the particle displacements are taken to calculate the mean square displacement variations in time and so the diffusion coefficients. To monitor the random motion of non-spherical particles a reference frame is used that the particle just have translational motion in it. We call it the body frame that is just like the particle rotates about the z-axis of the lab frame. Some statistical analysis, such as velocity autocorrelation function and histogram of displacements for the particle either in the lab or body frames, are available in the program. Program also calculates theoretical values of the diffusion coefficients for particles of some basic geometrical shapes; sphere, spheroid and cylinder, when other diffusion parameters like temperature and fluid viscosity coefficient can be adjusted. Program summary Program title: KOJA Catalogue identifier: AEHK_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEHK_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 48 021 No. of bytes in distributed program, including test data, etc.: 1 310 320 Distribution format: tar.gz Programming language: MatLab (MathWorks Inc.) version 7.6 or higher. Statistics Toolbox and Curve Fitting Toolbox required. Computer: Tested on windows and linux, but generally it would work on any computer running MatLab (MathWorks Inc.). There is a bug in windows 7, if the user is not the administrator sometimes the program was not able to overwrite some internal files. Operating system: Any supporting MatLab (MathWorks Inc.) v7.6 or higher. RAM: About eight times that of loaded data Classification: 12 Nature of problem: In many areas of physics, knowing diffusion coefficients is vital and gives useful information about the physical properties of diffusive particles and the environment. In many cases a diffusive particle is not a sphere and has rotation during its movements. In these cases information about a particle’s trajectory both in lab and body frame would be useful. Also some statistical analysis is needed to obtain more information about a particle’s motion. Solution method: This program tries to gather all required tools to analyse raw data from the Brownian motion of a diffusing particle. Ability to switch between different methods of calculation of mean square displacement to find diffusion coefficients depends on the correlations between data points. There are three methods in the program: time average, ensemble average and their combinations. A linear fit is

✩ This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/ science/journal/00104655). Corresponding author. E-mail address: [email protected] (M.A. Charsooghi).

*

0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.09.017

© 2010

Elsevier B.V. All rights reserved.

M.A. Charsooghi et al. / Computer Physics Communications 182 (2011) 400–408

401

done to measure Diffusion Coefficient (D), the weight and fraction of data points is controllable. Given physical properties of the system, the program can calculates D theoretically for some basic geometrical shapes; sphere, spheroid and cylinder. In the case of non-spherical particles if data of rotation is available, the code can calculate trajectory and diffusion also in body frame. There are more statistical tools available in the program, such as histogram and autocorrelation function to obtain more information e.g. relaxation time to ideal diffusion motion. Code uses log–log diagram of mean square displacement (MSD) to calculate the amount of deviation from normal diffusion to sub- or super-diffusion. Running time: It is dependent on the input data, but for typical data in the order of mega bytes, it would take tens of minutes. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Brownian motion of molecules of a fluid is an efficient mechanism for material transport inside the fluid by means of the diffusion processes. For example, complex functionality of living cells, depends on Brownian motion of their constituents. Detailed information about the environment can be gained by analysis of the particle’s trajectory [1–3]. The physical nature of interactions between a Brownian particle and its surrounding has been investigated in many research works [4–11]. Investigating the time evolution of the mean square displacement, gives us information about the diffusive behaviour of the particle. Different many-particle techniques have been widely used to study the diffusion processes. Dynamic light scattering [12–16], neutron scattering [17–19], and fluorescence correlation spectroscopy (FCS) [20] are from the most common ones. In recent years, optical methods have been developed to observe the motion of single particles [21,22]. These methods, called single particle tracking (SPT), and based on recording the trajectory of a moving particle. In SPT methods transport properties of the particle are derived by statistical analysis of the particle motion trajectories. Such method provides some information that are not available through measurements on many-particle systems. For example Saxton et al. measure trajectories of individual proteins or lipids in the plasma membrane of cells [22], Kusumi et al. observed movement, recruitment, and activation of single molecules in the plasma membrane in living cells [23], and Dahan et al. studied the diffusion dynamics of glycine receptors [24]. Diffusion behaviour of a spherical particles can be characterize by their translational diffusion coefficient D but for an asymmetric particle, in an isotropic homogeneous fluid, at least there are two different translational diffusion coefficients, D a for diffusion along the long or principal axis and D b along the minor axis [25–27], and a rotational diffusion coefficient, D θ [28,29]. Here we introduced a graphical user interface package which calculates the diffusion coefficients from experimental or simulated recorded particle trajectories. Also for asymmetric particles, the program calculates the motion trajectories in body frame (a reference frame which the particle just has linear translation in it) and computes translational and rotational diffusion coefficients in lab and body frame. Other statistical quantities such as velocity autocorrelation function (VAF) and histogram of displacement values of particle are also accessible in the program. In this program we tried to gather all needed tools for analysis of a Brownian motion of an (as)symmetric particle in a single user friendly interface. There are three methods in the program for these calculations: time average, ensemble average and their combinations. Ability to switch between different methods of calculation of mean square displacement (MSD) to find diffusion coefficients depends on the correlations between data points. Based on the concept of a diffusion process, the slope of a linear fit to the variations of the calculated MSDs of the particle versus time, provides the diffusion

coefficients. The program also is able to control the weight and the number of data points. For given physical properties of the system, the program can calculate D for particles of some basic geometrical shapes; sphere, spheroid and cylinder based on theoretical models by Perrin and Li et al. [28–30]. In the case of non-spherical particles if the data of particle orientation are available, the code can calculate trajectory and translational diffusion coefficients in body frame. In continuance, the physical basis of the calculations in the program will be discussed. Then some equations which are used to calculate theoretical values of D will be introduced. After that we will explain different parts of the programs in detail and show how user can use the program. At the end a summary of the program’s benefits and advantages are presented. 2. Data analysis The input data for the program corresponds to the positions and orientation of the particle recorded in an experiment or generated in a simulation process. One of the main challenges in experimental data is the number of data points. To have a good estimate of diffusion coefficients a large number of data points are required. Otherwise some tricks should be used to obtain trustable and precise results. If the data points are uncorrelated sub-ensemble can be defined and a combination of time average and average over the sub-ensembles will increase the precision of the results in comparison with the case when just ensemble average is taken over the whole data points. In the other hand using both methods separately and comparison of the results gives a chance to have some conclusion about correlations between the data points. In this program we have introduced functions for all of these situations which in this section will be explained in details. 2.1. Diffusion coefficients The Brownian diffusion coefficient, D of an isolated spherical particle is inversely proportional to the drag (or friction) coefficient γ via the Einstein relation,

D=

kB T

γ

.

(1)

Here kB is the Boltzmann constant and T is the temperature. In the case of Brownian motion of a particle, D can be measured by [31],

   2 r(t )2 − r(t ) = 2nDt ,

(2)

where r(t ) is equal to the displacement of the particle during a time interval of t, r(t ) − r(0), and n is the space dimension. In general it’s convenient to calculate D in each dimension separately, e.g. x, y or z, if exist. So Eq. (2) is simplified to,

   2 x(t )2 − x(t ) = 2D x t ,

(3)

402

M.A. Charsooghi et al. / Computer Physics Communications 182 (2011) 400–408

Fig. 1. The procedure of ensemble averaging over data, as an example the boxes with dashed borders shows the ensemble average procedure on particle displacement during the time interval t j − t 0 .

where D x is the diffusion constant along the x-axis and for a homogeneous space, x can be any of the coordinates. For an isolated spherical particles, diffusion coefficients in all directions should be the same, in spite of other cases such as ellipsoids. From now on, we will restricted ourselves to one dimension, but the results can be easily extended to higher dimensions (the program calculates diffusion coefficients for all available dimensions). 2.2. Mean square displacements, ensemble average Assume that for a single particle there are m trajectories of a one-dimensional Brownian motion obtained empirically, or from simulations. Each trajectory describes the particle motion in n discrete time steps. xi (t j ) is the x component of the position of the particle after jth time steps over the ith trajectory (1 < i < m and 0 < j < n). To obtain the ensemble average of the particle displacement during the time interval t j − 0, the program calculates xi (t j ) = xi (t j ) − xi (0) and (xi (t j ))2 and then calculate D according to Eq. (3) m 1 

m

2 xi (t j ) −

i =1



m 1 

m

where σ is the standard deviation of the samples and n is the number of samples. So if the number of trajectories are small, the number of samples which we average over to calculate each point of MSD diagram, would be small. Therefore the error on the final measured D would be large. In this case using time average method increases the number of samples in the averaging process. In mathematical words, in this case, displacement in time interval τ j , x(τ j ) − x(0), is similar to the displacement in the time interval between t j  and t j  + τ j , x(t j  + τ j ) − x(t j  ). We can consider τ j = j δt ( j = 0, 1, 2, . . . , n), which δt is the minimum time interval. In this case the time average of displacements in time intervals of τ j can be written as



n− j

  xi (t j  + τ j ) − xi (t j  ) xi (τ j ) = , n− j 

= 2Dt



i =1

for j = 0, 1, 2, . . . , n,

(4)

(6)

If there is any correlation between the data points, then the time averaging appeared in Eq. (6) will leads to incorrect results. Instead, the following relation should be used to eliminate the overlap between the motion intervals;

2

 xi (t j )

i = 1, 2, 3, . . . , m .

j =0

 n− j 

j   xi ((1 + j  )τ j ) − x( j  τ j ) xi (τ j ) = ,  n−j j  j  =0

i = 1, 2, 3, . . . , m , (7)

which averaging is over different trajectories. The second term in the left-hand side comes from the particle velocity drift or any flow in the system. If there isn’t any flow, it would be automatically zero. The ensemble averaging over the ensemble of trajectories is shown clearly in Fig. 1. At the end, mean square displacement (MSD), x2 (t ) − x(t )2 , as a function of time will be obtained. At each point of MSD diagram, error of calculated MSD, is standard deviation of averaging over different trajectories divide by square root of number of samples. Diffusion coefficient is deduced from the slope of a linear line which fitted to MSD diagram. 2.3. Mean square displacements, time average

σ

D=

m 

D j /m.

(8)

j =1

2.4. Mean square displacements, time-ensemble average

The standard deviation of the mean value is related to the number of samples by the relation

standard deviation of the mean = √ , n

where x is the integer part of x. Both averaging with and without overlapping is considered in the program. Figs. 2 and 3 show the difference between what we called correlated and uncorrelated time averaging method. Performing time average over the jth trajectory leads to a diffusion coefficient D j for this trajectory and finally the diffusion coefficient D, can be written as

(5)

This method is a combination of the two previous methods. The program first calculates (x(t + τ ) − x(t ))2 and (x(t + τ ) − x(t ))2  using time average over each trajectory, then take ensemble average of these two quantities over all different trajectories.

M.A. Charsooghi et al. / Computer Physics Communications 182 (2011) 400–408

403

Fig. 2. The procedure of time averaging over data. Arrows with solid, dash and dash-dot outlines show the time steps of t = 1δt, t = 2δt and t = 3δt respectively, which δt is the minimum time step in the system. A box with dash-dot border shows that the procedure of averaging is done only on separate data sets. A box with dashed boundary and gray background shows the area of overlapping between time intervals t 2 − t 0 and t 3 − t 1 for t = 2δt. In this method overlapping occurs for all time intervals, except t = 1δt.

Fig. 3. The procedure of time averaging over data without any overlapping. As an example rows with solid and dash outlines show the time steps of t = 1δt and t = 2δt respectively. A box with dash-dot border shows that the procedure of averaging is done only on separate data sets.

2.6. Autocorrelation function

2.5. Anomaly In a normal diffusion process, the mean squared displacement (MSD) of a particle is a linear function of time. Anomalous diffusion is used to describe a diffusion process with a non-linear dependence on time. Diffusion is often described by a power law,

MSD(t ) ∼ Dt α ,

(9)

where D is the diffusion coefficient, t is the elapsed time and α is known as diffusion anomaly. In a typical diffusion process, α = 1. If α > 1 (α < 1), the phenomenon is called super-diffusion (subdiffusion). As an example, super-diffusion can be observed as a result of active cellular transport processes. and sub-diffusion has been proposed as a measure of macromolecular crowding in the cytoplasm [32–34]. One can write log Dt α = log D + α log(t ) and to find α , it’s only enough to calculate the slope of the linear fit to variations of log(MSD) versus log(t ). Since the value of α has some errors, it is better to set a range on α to obtain the diffusion anomaly. For example instead of α = 1 for normal diffusion, it’s better to use |α − 1|  δ α , when δ α is an error in measured α .

The velocity autocorrelation function can be used to check if the particle motion is a simple diffusion or not. According to the theory of ideal diffusion, velocity autocorrelation function, defines as [35],





v (t ) v (0) = 2D δ(t ),

(10)

where δ(t ) is the Dirac delta function and v (0) is velocity of the particle at t = 0. In the program, velocity is simply defined from data of position and time using equation

v (i ) =

x(i + 1) − x(i ) t ( i + 1) − t ( i )

,

i = 0, 1 , 2 , . . . , n − 1 ,

(11)

where i stands for the ith time step of motion. When the time average is taking over a time interval that is much larger than inter molecular collision time (i.e. < 10−12 s), Eq. (10) is well satisfied and quite clear. Large number of collisions during the measurement time cause the velocities of the particle in successive time intervals become uncorrelated [36]. For a particle with simple diffusion, the velocity autocorrelation function is computed using relation





A (t ) = v (t + τ ). v (τ ) ,

(12)

404

M.A. Charsooghi et al. / Computer Physics Communications 182 (2011) 400–408

Fig. 6. Coordinate transformation for rotation of a spheroid. The primed coordinate fixes to spheroid axes and rotates with rotating of a spheroid. θn is the angle between X directions of lab and body frames, in nth step.

2.8. Trajectory and diffusion in the particle body frame Fig. 4. The velocity autocorrelation function versus time for random motion of a 3μ spherical polystyrene bead inside water at room temperature; experiments – solid line, exponential fit – dashed line.

For non-spherical particles, e.g. a spheroid with long axis of length 2a and short axes of length 2b, translational diffusion is anisotropic. We are considering that the particle only has a planar rotation and defining a new reference frame ( X  , Y  , Z  ) that just rotates with the particle about the Z -axis of the lab frame when Z = Z  . We call this reference frame as the body frame. In the body frame, the particle long (short) axis is along the X  (Y  ) direction and the particle diffusion along its long (short) axis is



specified by diffusion coefficient, D t = k B T /γt (D t⊥ = k B T /γt⊥ )

Fig. 5. The histogram of displacements for random motion of a 3μ spherical polystyrene bead inside water at room temperature; experiments – dashed line, exponential fit – solid line.

where A (t ) is velocity autocorrelation function. Eq. (12) should have a non-zero value for t = 0. By fitting an exponential function, −t

e β to A (t ), we can find β which specifies a time constant for a simple diffusion. When t < β the diffusion is not simple. Typically β has small values but if any extra forces beside of the thermal collisions influences the random motion of the molecules, it would be grow. In Fig. 4 a typical autocorrelation function for a random walker is shown.

(see Fig. 6), when γt and γt⊥ are viscous friction coefficients parallel and perpendicular to the main axis, respectively. The particle rotation can be specified by θ , the smaller angle between X and X  axes and the rotational diffusion coefficient is defined by D r = k B T /γr where γr is rotational viscous friction coefficient. Perrin analytically calculated the viscous friction coefficients, for spheroid in three dimensions [28,29]. Displacement of the particle in nth step of its random motion in x direction is dx(tn ) = x(tn+1 ) − x(tn ) and in y direction is dy (tn ) = y (tn+1 ) − y (tn ) (Fig. 6). To find the trajectory of the particle in the body frame, the displacements dx(tn ) = x(tn+1 ) − x(tn ) and dy (tn ) = y (tn+1 ) − y (tn ) and also changes in rotation angle by dθ(tn ) = θ(tn+1 ) − θ(tn ) are calculated for each time step. Displacements in body frame (with primed symbols) are related to lab frame displacements via



dxn



dyn

cos θn sin θn − sin θn cos θn



dxn dyn

,

(14)

where θn = (θ(tn ) + θ(tn+1 ))/2 [37]. Then in the body frame one can write the position for the time t j as:

x (t j ) =

j 

dxi ,

j = 1, 2, 3, . . . , n ,

i =1

y  (t j ) =

2.7. Histogram

=

j 

dy i ,

j = 1, 2, 3, . . . , n .

(15)

i =1

This program also calculates histogram of displacements with a Gaussian fit to it. Width of the Gaussian distribution which is fitted to histogram of displacement is related to diffusion coefficient via 2

σ = 2Dt , 2

(13)

where σ is the variance of the distribution and, D is the diffusion coefficient. Fig. 5 shows a typical diagram of displacements histogram which is produced with the program.

2.9. Theoretical values of diffusion coefficients To find some estimation of the diffusing particle shape or to test how the calculation presumptions matches with the theoretical diffusion models, the program also calculates the theoretical values of diffusion coefficients for particles with one of the three basic geometrical shapes; sphere, spheroid and cylinder. The diffusion coefficient, D, and the viscous friction coefficient are related

M.A. Charsooghi et al. / Computer Physics Communications 182 (2011) 400–408

according to Einstein’s relation (Eq. (1)). For a spherical shaped particle, γ is equal to [31]

γr =

γ = 6πηa

γr⊥ =

(16)

where η is the viscosity of surrounding media and a is the radius of the sphere. For a spheroid with semi-axis length of a and b = c,

translational viscous friction coefficients along, γt , and perpendic⊥ ular, γt , to the axis a are [28,29],

γt = 16πη

a2 − b 2

(17)

, − b2 )s − 2a a2 − b 2 γt⊥ = 32πη 2 . (2a − 3b2 )s + 2a (2a2

(18)

Typically there are two different spheroids, prolated spheroid and oblated spheroid. In prolated spheroid a > b and s is

2

s= √ log a2 − b 2

a+

√

a2 − b 2 b

,

(19)

and for oblated spheroid a < b and s is

2

s= √ arctan b 2 − a2

b

,

(2a2

(20)

and for rotation about the mentioned axes the rotational friction coefficients are,

− b2 )s − 2a

32πη (a2 − b2 )b2 2a − b2 s

3

(21)

,

(22)

.

For cylindrical particles with diameter d and length L and so the axial ratio of p = L /d the translational viscous friction coefficients parallel and perpendicular to the principle axis and rotational coefficients around these axes are [30]

γt = γt⊥ = γr = γr⊥ =

2πη L

(longwise),



ln p + δt 4πη L

(23)

(sidewise),

ln p + δt⊥

(24)



3.84πη L 3 (1 + δr ) p2

πη L 3 3(ln p + δr⊥ )

√

b 2 − a2

a4 − b 4

32πη 3

405

(about long axis),

(about short axis),

(25) (26)



where δt , δt⊥ , δr , and δr⊥ are the end correction coefficients which are functions of axial ratio, p. For cylinders with p > 1 but not for p = ∞ limit we will have [30,38]:

δt = −0.207 +

Fig. 7. Flowchart of the program structure.

0.980 p

−

0.133 p2

,

(27)

406

M.A. Charsooghi et al. / Computer Physics Communications 182 (2011) 400–408

δt⊥ = 0.839 +

δr =

0.677 p

−

0.185 p 0.183

δr⊥ = −0.662 +

p2

+

0.233 p2

(28)

,

(29)

,

0.917 p

−

0.05 p2

(30)

.

3. Using the program To work with program open KOJA.m in MATLAB 7.6 or higher and run the m-file. A brief review of the program structure, as a flowchart, is shown in Fig. 7. Different parts of the program are introduced in the following sub-sections. Each input file should corresponds to Brownian motion of a single particle for a specified time duration and time steps when at each time step the x, y and z component of the particle position and the particle orientation respect to x-axis in the laboratory frame, θ , are recorded. The input data file is in the form of a matrix when the elements on each row represents the particle position and orientation in a specified time step and the columns corresponds to the time, x, y, z position coordinates and the particle orientation θ . This means that the maximum number of the columns are five and number of the rows are determined by the time duration of the motion and length of time steps. The program also recognizes the user defined headers for input files. The typical structure for input file of a nTable 1 Typical structure of an input file. Each row corresponds to a time step and columns show particle’s position and orientation in associated time. Time

Coordinates

t

x

y

z

θ

t1 t2 t3

x1 x2 x3

y1 y2 y3

z1 z2 z3

. . .

. . .

. . .

. . .

θ1 θ2 θ3 . . .

time step Brownian motion is shown in Table 1. In each run, the structure of the input files should be the same. 3.1. Windows There is a Main Window with four sub-windows which user deals with. 3.1.1. Main window After running the program, a window which we call it “Main window” appears (Fig. 8). It contains some major parts: Axes Area, Current Frame, Durrent Data, Fit Options panels. Axes Area – Each time user select one of sub-menus, toolbar icons or any item of control panels (the panels left to the Axes Area), the related physical quantity is computed and associated plot will appear in this part. User can save the resulting figure using file menu or toolbar icon, in png, bmp, eps, tiff, jpg and other graphical formats. Current Frame – User can switch between Body Frame and Lab Frame, using this panel. Current data – In this part user can switch between the motion Dimensions and the Data sets. After selection, the data which has shown in the plotting area, will be automatically refreshed. Fit Options – In this part user can change Fit Options. User can choose a portion of the data starting from the beginning of the data points and running the fitting procedure based on the linear fit. It is also possible to set the Weight of each data point in the this procedure. Like in the MATLAB fitting procedure, the weight of each data point is assigned according to the error value at that point. The fitting weight, w, is written as,

Orientation

w=

1

σ2

(31)

where σ is the error on each data point. Program uses exponential and Gaussian fit for velocity autocorrelation function and histogram, respectively.

Fig. 8. A screen shot of a main window of the program. User has access to all applications and functions through menu-bars and toolbars.

M.A. Charsooghi et al. / Computer Physics Communications 182 (2011) 400–408

Fig. 9. When user click α sign from toolbar icons or click Calculate Theory from Diffusion menu, the theory sub-window will be open. User can calculate translational and rotational diffusion for three basic geometrical shape; sphere, spheroid and cylinder.

3.1.2. Theory sub-window When user click α sign from toolbar icons or click Calculate Theory from Diffusion menu, the theory sub-window will be open. Fig. 9 shows a snapshot of a theory sub-window. User can calculate translational and rotational diffusion for some basic geometrical shapes; sphere, prolate spheroid, oblate spheroid and cylinder. User can set sizes, units, temperature, and viscosity of medium. 3.1.3. Results sub-window When user click Results from View menu, the results subwindow will be opened. User can follow, save, edit and clear the results of calculations there. Especially for diffusion calculations, program exports diffusion and velocity of drift with their amount of error to this window. 3.1.4. Description sub-window When user click Description from View menu, the description sub-window will be opened. User can write any extra information about data, such as conditions of experiment or simulation, there. Also user can edit, save and clear it and load any saved description. 3.1.5. Log sub-window When user click Log from View menu, the log sub-window will be opened. User can follow, save, edit and clear its activity through this window. 3.2. Menus 3.2.1. File menu User can load raw data, open and save its projects and figures and also export axis data or print it using this menu. 3.2.2. Diffusion menu Time, Ensemble, Time-Ensemble Average – Each of these submenus calculates diffusion coefficient using the associated method. These methods are explained completely in Section 2. MSD versus time, with a linear fit will be shown in the figure panel. The value of the coefficients are written in Results panel. Anomaly – To calculate non-linearity dependence of diffusion on time, user should use this sub-menu. Before that user should

407

Fig. 10. Diffusion statistics. When user use time average method to calculate diffusion coefficient. There would be diffusion coefficients as many as the number of data sets. Use this sub-menu to plot these coefficients together. Each point corresponds to one diffusion coefficient. Dashed line show the final diffusion coefficient.

plot MSD versus time in the figure panel using one of the methods of calculating diffusion coefficients (Time, Ensemble, TimeEnsemble average). Then the log–log plot of MSD and time is plotted and the kind and amount of anomaly will be presented in Results panel. Diffusion Statistics – If diffusion coefficients compute using time average method, there will be a diffusion coefficient for each data set. To comparison these coefficients using this sub-menu will plot all the values in one figure. The number of points is equal to number of data sets (Fig. 10). If user check the Show Fit box, the final value of diffusion will be add as a straight line. Calculate Theory – When user click α sign from toolbar icons or click Calculate Theory from Diffusion menu, the theory subwindow will be open. For more information refer to part 3.1.2. 3.2.3. Statistics menu This menu will calculate and show Histogram of displacements and velocity autocorrelation function. User can fit Gaussian function to histogram and exponential function to autocorrelation diagram. Fit parameters is shown in Results panel. 3.2.4. View menu This menu contains some visual commands to change the face of the program’s panels. User can see results, description and log sub-windows. Also user can plot trajectory of particle motion versus time or in space. 3.2.5. Plot menu User can change the view of the plotting area switching between simple, log–log, semilogx and semilogy. There are options for showing errors and grids both major and minor on the plot. Using Edit Axes options, the figure is shown in the new window with many extra plot tools such as colors, labels, thicknesses, markers and etc. User can also clear the axes area using this menu. 3.2.6. Help menu Product Help – A user manual for this package in pdf format. About – About the version of this product and a collaboration team which produce the software. User can save all calculated quantities as a mat file. All calculated data is loaded any time user opens its saved project. As diffu-

408

M.A. Charsooghi et al. / Computer Physics Communications 182 (2011) 400–408

sion parameters are important parts of this program, we mention them separately. The value of diffusion, its error and its intercept, velocity and its error and the evolution of MSD can be saved. In the case of calculation of diffusion coefficients using time average method, beside the parameters of final case, parameters for each trajectory are also available and can be saved. These situations are established for both lab and body frame. For more information about using the program, user can refer to the user manual which is available with the program.

Acknowledgements We are grateful to Farshid Mohammad-Rafiee, Sarah Mohammadinejad, Laleh Mollazadeh and Faegheh Hajizadeh for helpful discussions. Appendix A. Supplementary material Supplementary material related to this article can be found online at doi:10.1016/j.cpc.2010.09.017.

4. Program features References In the following some useful features of the program are highlighted

• The program can call all data files recorded in txt, dat, and

•

•

•

•

•

•

ascii formats. The program also benefits from a user friendly menu bar and sub-menus to manipulate different program windows. Since sometimes specially in the simulation processes, a large amount of data are generated, to optimize memory usage and prevent any out of memory error, the program saves all variable values on hard disk. For large amount of data (few hundred mega bytes), calculations take long time and user should be patient. It’s a graphical user interface (GUI) program, so helps user to have everything in its hands. This prevents producing too many output figures as the results of different calculations on the loaded data. To observe different outputs it is just enough to switch to different applications on the menu bar like “Trajectory” for plotting the displacement in x, y or z direction versus time or “MSD” to find the mean square displacement in x, y or z directions. Whole of a project can be saved so user can save all calculated parameters and also the next time when the project is opened again, user will have access to all the data and results, which calculated in last session. In the first running user should determine the physical quantity of each column of data, but after that the program will memorize until structural changes appear in data files. Program consider user’s first choices as defaults but user can changes them in any run. A report of the output results and a log file in a txt format also are available. User can take a report and also log file of his/her activities in the program. Also user can save, clean and upload them. This program is not only applicable for diffusion of colloids or microorganisms in fluids, but also it can analyse any kind of data of diffusion processes obtained from experiments or computer simulations.

[1] F.C. MacKintosh, C.F. Schmidt, Current Opinion Colloid Interface Sci. 4 (1999) 300. [2] C. Tischer, et al., Appl. Phys. Lett. 79 (2001) 3878. [3] S. Jeney, E.H.K. Stelzer, H. Grubmüller, E.-L. Florin, Chem. Phys. Chem. 5 (2004) 1150. [4] A. Einstein, Ann. Phys. 14 (2005) 182. [5] E.J. Hinch, J. Fluid Mech. Digital Archive 72 (1975) 499. [6] A. Rahman, Phys. Rev. 136 (1964) A405. [7] K. Ohbayashi, T. Kohno, H. Utiyama, Phys. Rev. A 27 (1983) 2632. [8] D.A. Weitz, D.J. Pine, P.N. Pusey, R.J.A. Tough, Phys. Rev. Lett. 63 (1989) 1747. [9] Y.W. Kim, J.E. Matta, Phys. Rev. Lett. 31 (1973) 208. [10] A. Meller, et al., Biophysical Journal 74 (1998) 1541. [11] E.J.G. Peterman, M.A. van Dijk, L.C. Kapitein, C.F. Schmidt, Rev. Scientific Instruments 74 (2003) 3246. [12] B.J. Berne, R. Pecora, Dynamic Light Scattering: With Applications to Chemistry, Biology, and Physics, Dover Publications, 2000. [13] S. Roldán-Vargas, et al., Phys. Rev. E 80 (2009) 021403. ˇ [14] A. Mertelj, L. Cmok, M. Copiˇ c, Phys. Rev. E 79 (2009) 041402. [15] P. Holmqvist, G. Nägele, Phys. Rev. Lett. 104 (2010) 058301. [16] E. Andablo-Reyes, P. Díaz-Leyva, J.L. Arauz-Lara, Phys. Rev. Lett. 94 (2005) 106001. [17] J.S. Higgins, H.C. Benoît, Polymers and Neutron Scattering, Oxford Series on Neutron Scattering in Condensed Matter, vol. 8, Oxford University Press, USA, 1997. [18] T. Spehr, et al., Phys. Rev. E 79 (2009) 031404. [19] E. Mamontov, Y.A. Kumzerov, S.B. Vakhrushev, Phys. Rev. E 72 (2005) 051502. [20] K. Bacia, P. Schwille, Methods 29 (2003) 74. [21] H. Geerts, et al., Biophysical Journal 52 (1987) 775. [22] M.J. Saxton, K. Jacobson, Annual Rev. Biophys. Biomolecular Struct. 26 (1997) 373. [23] A. Kusumi, et al., Annual Rev. Biophys. Biomolecular Struct. 34 (2005) 351. [24] M. Dahan, et al., Science 302 (2003) 442. [25] P.G. Saffman, M. Delbrück, Proc. Nat. Acad. Sci. USA 72 (1975) 3111. [26] Y. Gambin, et al., Proc. Nat. Acad. Sci. USA 103 (2006) 2098. [27] A.J. Levine, T.B. Liverpool, F.C. MacKintosh, Phys. Rev. E 69 (2004) 021503. [28] F. Perrin, J. Phys. Radium 5 (1934) 497. [29] F. Perrin, J. Phys. Radium 7 (1936) 1. [30] G. Li, J.X. Tang, Phys. Rev. E 69 (2004) 061921. [31] P. Nelson, Biological Physics: Energy, Information, Life, W.H. Freeman, 2003. [32] A. Caspi, R. Granek, M. Elbaum, Phys. Rev. E 66 (2002) 011916. [33] M. Weiss, M. Elsner, F. Kartberg, T. Nilsson, Biophysical Journal 87 (2004) 3518. [34] J. Bouchaud, Phys. Reports 195 (1990) 127. [35] R. Zwanzig, Annual Rev. Phys. Chem. 16 (1965) 67. [36] S. Chandrasekhar, Rev. Modern Phys. 15 (1943) 1. [37] Y. Han, et al., Science 314 (2006) 626. [38] J.G. de la Torre, V.A. Bloomfield, Quart. Rev. Biophys. 14 (1981) 81.