GLE main file

Compute various physical quantities (MSD, IncACF, FPT) for a spherical particle freely diffusing in a Generalize Rouse model. The particle's motion is described using a GLE.

Warning: no path data is saved.

Christel Hohenegger

Last updated: 06/01/2016

Initialization
Generate data using SimManager
Plots
Test the goodness of the quadratic fit

Initialization

Parameters:

$\alpha$ is the Rouse exponent
radius is the radius of the particle in microns
width is the width of the layer in microns
number of time steps and time steps for the ACF are given in milliseconds

Outputs:

MSDMat: MSD array (time, alpha, radius, kernels)
VelACFMat: Inc ACF array (lag, alpha, radius, kernels)
stdFPTMat: standard deviation of the first passage time (used for error bars plotting)
FPTMat: FPT array (width, time, alpha, radius, kernels)
Survivors: Relative number of survivors (width, time, alpha, radius, kernels)

clearvars
close all

alpha = [2];
numKernels = [4 16];
radius = [1 0.5 2];
numTsteps = 2^13;
numTACF = 2^3;
width = (0.1:0.05:2);

% Allow the user to decide to plot the mean square displacement (MSD), the
% increment autocorelation (ACF) or the mean first passage time (FPT)
promptMSD = input('Do you want to plot the MSD [y]/n? ','s');
if isempty(promptMSD)
    promptMSD = 'y';
end
promptIncACF = input('Do you want to plot the increment ACF [y]/n? ','s');
if isempty(promptIncACF)
    promptIncACF = 'y';
end
promptMFPT = input('Do you want to plot the mean FPT [y]/n? ','s');
if isempty(promptMFPT)
    promptMFPT = 'y';
end
fprintf('\n')

MSDMat = NaN(numTsteps,length(alpha),length(radius),length(numKernels));
VelACFMat = NaN(numTACF-1,length(alpha),length(radius),length(numKernels));
meanFPTMat = NaN(length(width),length(alpha),length(radius),...
    length(numKernels));
stdFPTMat = NaN(length(width),length(alpha),length(radius),...
    length(numKernels));
FPTMat = NaN(length(width),numTsteps+1,length(alpha),length(radius),...
    length(numKernels));
Survivors = NaN(length(width),numTsteps+1,length(alpha),length(radius),...
    length(numKernels));

Generate data using SimManager

Order of operations:

define a general sim with the given parameters
generate the other physical parameters, override them is necessarily
calculate the residue
generate the covariance matrix
generate the paths
find the quantity of interest

for iindex = 1:length(alpha)
    alphaVal = alpha(iindex);
    fprintf('Working on i: %d.\n',iindex)
    for jindex = 1:length(radius)
        radiusVal = radius(jindex);
        fprintf('Workind on j: %d.\n',jindex)
        for kindex = 1:length(numKernels)
            fprintf('Working on k: %d.\n',kindex)
            numK = numKernels(kindex);
            if (strcmp(promptMSD,'y')||strcmp(promptIncACF,'y'))
                sim = SimManager(radiusVal,numK,alphaVal,numTsteps);
                sim.GenerateParameters;
                sim.parameters.dt = 0.01;
                sim.parameters.time = (0:sim.parameters.numTsteps)...
                    *sim.parameters.dt;
                sim.ResidueCalc;
                sim.GenerateCovMat;
                sim.GeneratePaths;
                if strcmp(promptMSD,'y')
                    sim.FindMSD;
                    MSDMat(:,iindex,jindex,kindex) = sim.dataMSD;
                else
                    sim.FindVelACF(numTACF);
                    VelACFMat(:,iindex,jindex,kindex) = sim.dataVelACF;
                end
            end
            if strcmp(promptMFPT,'y')
                numTsteps = 2^13;
                sim = SimManager(radiusVal,numK,alphaVal,numTsteps);
                sim.GenerateParameters;
                sim.parameters.numPaths = 3e3;
                sim.ResidueCalc;
                sim.GenerateCovMat;
                sim.GeneratePaths;
                sim.FindMFPT(width);
                FPTMat(:,:,iindex,jindex,kindex) = sim.dataFPT;
                Survivors(:,:,iindex,jindex,kindex) = sim.dataSurvivors;
                meanFPTMat(:,iindex,jindex,kindex) = sim.dataMFPT;
                stdFPTMat(:,iindex,jindex,kindex) = sim.dataSMFPT;
            end
        end
    end
end

Plots

Legends might need to be changed. List of figures

MSD for different number of kernels (default radius and $\alpha$ )
MSD for different $\alpha$ and number of kernels (default radius)
MSD for different radius and number of kernels (default $\alpha$ )
IncACF for different $\alpha$ and number of kernels (default radius)
MFPT for different number of kernels (default $\alpha$ and radius) rescaled to minutes
MFPT for different number of kernels (default $\alpha$ and radius) rescaled to minutes with 95\% confidence intervals (commented out)
MFPT for different number of kernels and radius (default $\alpha$ ) rescaled to minutes with 95\% confidence intervals (commented out)

if strcmp(promptMSD,'y')
    time = sim.parameters.time;
    figure(1)
    loglog(time(2:end),squeeze(MSDMat(:,1,1,:)),'LineWidth',1.5)
    xlabel('Time [ms]','FontSize',12)
    ylabel('MSD [micron^2]','FontSize',12)
    set(gcf,'color','w')
    set(gca,'FontSize',12)
    legend('N = 1','N = 5', 'N = 20', 'N = 40','Location','NorthWest')
    xlim([0 max(time)])
    figure(2)
    loglog(time(2:end),squeeze(MSDMat(:,:,1,1)),'LineWidth',1.5)
    hold all
    loglog(time(2:end),squeeze(MSDMat(:,:,1,4)),'LineWidth',1.5)
    xlabel('Time [ms]','FontSize',12)
    ylabel('MSD [micron^2]','FontSize',12)
    set(gcf,'color','w')
    set(gca,'FontSize',12)
    legend( '\alpha = 2, N = 1', '\alpha = 4, N = 1', ...
        '\alpha = 2, N = 40', '\alpha = 4, N = 40',...
        'Location','NorthWest')
    xlim([0 max(time)])
    figure(3)
    loglog(time(2:end),squeeze(MSDMat(:,1,:,1)).*repmat(radius,numTsteps,...
        1),'LineWidth',1.5)
%    loglog(time(2:end),squeeze(MSDMat(:,1,:,1)),'LineWidth',1.5)
    hold all
     loglog(time(2:end),squeeze(MSDMat(:,1,:,4)).*repmat(radius,numTsteps,...
         1),'LineWidth',1.5)
%    loglog(time(2:end),squeeze(MSDMat(:,1,:,4)),'LineWidth',1.5)
    xlabel('Time [ms]','FontSize',12)
    ylabel('MSD [micron^2]','FontSize',12)
    set(gcf,'color','w')
    set(gca,'FontSize',12)
    legend('r = 1, N = 1', 'r = 0.5, N = 1', 'r = 2, N = 1', ...
        'r = 1, N = 40','r = 0.5, N = 40', 'r = 2, N = 40', ...
        'Location','NorthWest')
    xlim([0 max(time)])
end

if strcmp(promptIncACF,'y')
    figure(4)
    timeACF = (0:(numTACF-2));
    plot(timeACF,squeeze(VelACFMat(:,1,1,:)),'-o',...
        'MarkerSize',8,'LineWidth',1.5,'MarkerFaceColor','auto')
    hold
    plot(timeACF,squeeze(VelACFMat(:,2,1,:)),'-o',...
        'MarkerSize',8,'LineWidth',1.5,'MarkerFaceColor','auto')
    xlabel('Time Lag [ms]','FontSize',12)
    ylabel('Increment ACF','FontSize',12)
    set(gcf,'color','w')
    set(gca,'FontSize',12)
    legend('\alpha = 2, N = 1','\alpha = 2, N = 5', '\alpha = 2, N = 20',...
        '\alpha = 2, N = 40','\alpha = 4, N = 1','\alpha = 4, N = 5', ...
        '\alpha = 4, N = 20','\alpha = 4, N = 40','Location','NorthEast')
end

if strcmp(promptMFPT,'y')
    errorMat = 2.8*stdFPTMat/sqrt(sim.parameters.numPaths);
    errorMat = errorMat*1e-3/60;
    meanFPTMatMin = meanFPTMat*1e-3/60;
    figure(5)
    plot(width,squeeze(meanFPTMatMin(:,1,1,:)),'-o')
    errorbar(repmat(width',1,length(numKernels)),...,
        squeeze(meanFPTMatMin(:,1,1,:)),squeeze(errorMat(:,1,1,:)))
    xlabel('Width [micron]','FontSize',12)
    ylabel('Mean First Passage Time [minutes]','FontSize',12)
    set(gcf,'color','w')
    set(gca,'FontSize',12)
    xlim([0.05 2.05])
    legend('N = 1','N = 2', 'N = 4', 'N = 8','N = 16','Location','NorthWest')
%{
    figure(6)
    plot(width,squeeze(meanFPTMat(:,2,1,:)),'-o')
    errorbar(repmat(width',1,4),squeeze(meanFPTMat(:,2,1,:)),...
        squeeze(errorMat(:,2,1,:)))
    xlabel('Width [micron]','FontSize',12)
    ylabel('Mean First Passage Time [minutes]','FontSize',12)
    set(gcf,'color','w')
    set(gca,'FontSize',12)
    legend('N = 1','N = 5', 'N = 20', 'N = 40','Location','NorthWest')
    figure(7)
    plot(width,meanFPTMat(:,1,:,1),'-o')
    errorbar(repmat(width,1,4),squeeze(meanFPTMat(:,1,:,1)),...
        squeeze(errorMat(:,1,:,1)))
    hold all
    plot(width,meanFPTMat(:,1,:,4),'-o')
    errorbar(repmat(widcloseth,1,4),squeeze(meanFPTMat(:,1,:,4)),...
        squeeze(errorMat(:,1,:,4)))
    xlabel('Width [micron]','FontSize',12)
    ylabel('Mean First Passage Time [minutes]','FontSize',12)
    set(gcf,'color','w')
    set(gca,'FontSize',12)
    legend('r = 0.5, N = 1', 'r = 1, N = 1', 'r = 2, N = 1', ...
        'r = 0.5, N = 40','r = 1, N = 40', 'r = 2, N = 40', ...
        'Location','NorthEast')
%}
end

Test the goodness of the quadratic fit

Test the goodness of the quadratic fit of the MFPT as a function of the width on log-log scale (commented out)

%{
for iindex = 1:length(alpha)
    fprintf('Working on i: %d.\n',iindex)
    for jindex = 1:length(radius)
        fprintf('Workind on j: %d.\n',jindex)
        for kindex = 1:length(numKernels)
            fprintf('Working on k: %d.\n',kindex)
            [logFit,gof]=fit(log(width)',...
            squeeze(log(meanFPTMatMin(:,iindex,jindex,kindex))),'poly1');
            gof.rmse
            logFitVal(:,iindex,jindex,kindex) = logFit(log(width));
            polyFitVal = exp(logFitVal);
            polyFitCoeff(1,iindex,jindex,kindex) = exp(logFit.p2);
            polyFitCoeff(2,iindex,jindex,kindex) = logFit.p1;
        end
    end
end
%}

GLE main file

Contents

Initialization

Generate data using SimManager

Plots

Test the goodness of the quadratic fit