Shortest paths PBC calculation

Shortest paths calculation in periodic CGSD simulation

Introduction

This example shows how to compute the shortest path between a bead in a CGSD simulation and its periodic image using the NetworCh library. The scripts and data files corresponding to this tutorial can be found in directory <NetworChDir>/tests/shortest_path_pbc/.

Loading and replicating the structure

The chain structure corresponding to a CGSD simulation is outputted in two files: a beads position file (e.g. <simulationName>.cn) which contains the id and the spatial coordinates of every beads in the simulation, and a bond file (e.g. <simulationName>.bnd) which contains the connectivity information between the beads. The structure can be loaded in Matlab using function read_cgsd by specifying the name of the simulation cn and bnd files (without their extension). For instance, files 0_percent_strain_x500.cn and 0_percent_strain_x500.bnd can be loaded with:

% Load CGSD structure
filename = '0_percent_strain_x500';
[cn,bnds,vol] = read_cgsd(filename);

By default, a cross-link is represented by a cross-link bond of zero length between two beads. To simplify the structure and reduce the computation time, one can merge the cross-linked beads into a single bead by using the merge_cross_links function:

% Merge cross-linked nodes
[cn,bnds] = merge_cross_links(cn,bnds);

The physical length of each bond can be calculated using function all_segments_length:

% Calculate bonds physical lengths
bonds_length = all_segments_length(cn,bnds,vol);

Since we are interseted in calculated the shortest path from a bead to its periodic image, we need to replicate our structure by periodicity. For instance, if we wish to calculate the shortest path between a bead and its replica along the y-direction, we can use function replicate_volume as follows:

ndir = 2; % replicate in the y-direction
[cn1,bnds1,vol1] = replicate_volume(cn,bnds,vol,ndir);

where beads position array cn1 will contain both the position of the beads in the primary volume vol in cn1(1:size(cn,1),:) and their replica in the y-direction in cn1(size(cn,1)+1:end,:) such that size(cn1,1) = 2*size(cn,1).

bonds_weight = bnds1(:,4); %use chain length as bond weight

% or

bonds_weight = bnds1(:,end); %use physical distance as bond weight

Computing the periodic image shortest path distribution

n = size(cn,1);

np = ne-ns+1;
sp = zeros(np,1);

conn1 = generate_connectivity(size(cn1,1),bnds1);

A1 = adjacency(cn1,bnds1,bonds_weight);
[~,comp1] = find_components(conn1);

A1S = sparse(A1);

for i=1:np
    if type == 2
        j = nid(i);
    else
        j = i+ns-1;
    end
    if bio_graph
        [d,~,~] = graphshortestpath(A1S,j,j+n);
        sp(i) = d;
    else
        [d,~] = shortest_path(conn1,A1,j,j+n,comp1);
        sp(i) = d;
    end
end

% Plot shortest paths distribution (only for cross-linked nodes)
nid = cn(nid,end-1);
output_data = [nid,sp,cn(:,end)];
sp = output_data(output_data(:,3)==1,2);

figure(1);
hist(sp,20);
xlabel('shortest path length');