why the variables in the parfor loop are not sliced or temporary variables in the Monte Carlo simulation?

Dearest Matlab users,

I am dealing for the first time with the parfor command. I have a Monte Carlo simulation code to simulation the materials evolution, it runs a lot of time so I would like to run it faster through multithreading. First I built a matrix contains 40000 sites, for every monte carlo step, 40000 attemps should be gone through.The following is the first my parfor-loop code. when I run the whole programme, error shows Error using ==> parallel funtion in line 598 Error in ==> parallel_function>make_general_channel/channel_general at 894 Index exceeds matrix dimensions.and If I change D2, it says error: The temporary variable D2 in a parfor is uninitialized.or D2 is indexed but not sliced. what should I do to change B2 or D2 to a matrix that can be valid in parfor loop?? Thanks so much!

N=200; Q=50; N_mcs =100; PT=40; %total (A+B+C), PT1=20; % PA<x<=PT1 belongs crystal B, PT1<x<=PT belongs to crystal A PA=0; % amorphous C %%%%%%%%%%%initialization%%%%%%%%%%%%%%%%%%%%%%%%% rand('state',110) A1=round(rand(1,N*N)*(PT-1)+1); %A1 or B1 is the matrix of species. A1=reshape(A1,N,N); A1=A1'; B1(3:2+N,3:2+N)= A1(1:N,1:N); B1(3:2+N,1:2)= A1(1:N,N-1:N); B1(3:2+N,N+3:N+4)= A1(1:N,1:2); B1(1:2,3:2+N)= A1(N-1:N,1:N); B1(N+3:N+4,3:2+N)= A1(1:2,1:N); B1(1:2,1:2)= A1(N-1:N,N-1:N); B1(1:2,N+3:N+4)= A1(N-1:N,1:2); B1(N+3:N+4,1:2)= A1(1:2, N-2+1:N); B1(N+3:N+4,N+3:N+4 )= A1(1:2,1:2); PercentageC = sum(sum(A1<=PA))/(N*N),%%check the percentage of amorphous material. PercentageB = sum(sum(PA<A1&A1<=PT1))/(N*N),%%check the percentage of crystal B material. PercentageA = sum(sum(PT1<A1&A1<=PT))/(N*N),%%check the percentage of crystal A material. [I,J] = find(A1>PA&A1<=PT1); %find the sites of the crystal B. NN=length(I); A2=zeros(N,N);% the sites the amorphous material occupy was assigned to 0. AA2=round(rand(1,NN)*(Q-1)+1); %A2 or B2 is the matrix of spin number(orientation). for i=1:NN A2(I(i),J(i))=AA2(i); end [K,L] = find(PT1<A1&A1<=PT); %find the sites of the crystal A. MM=length(K); AA3=round(rand(1,MM)*(Q-1)+51); %A2 or B2 is the matrix of spin number(orientation). for j=1:MM A2(K(j),L(j))=AA3(j); end B2(3:2+N,3:2+N)= A2(1:N,1:N); B2(3:2+N,1:2)= A2(1:N,N-1:N); B2(3:2+N,N+3:N+4)= A2(1:N,1:2); B2(1:2,3:2+N)= A2(N-1:N,1:N); B2(N+3:N+4,3:2+N)= A2(1:2,1:N); B2(1:2,1:2)= A2(N-1:N,N-1:N); B2(1:2,N+3:N+4)= A2(N-1:N,1:2); B2(N+3:N+4,1:2)= A2(1:2,N-2+1:N); B2(N+3:N+4,N+3:N+4 )= A2(1:2,1:2); %rand('state',100) D2=cell(204,1); for i=1:204 D2{i}=B2(i,:) end D1=cell(204,1); for i=1:204 D1{i}=B1(i,:) end m=round(rand(N)*(N-1)+3); n=round(rand(N)*(N-1)+3); s=reshape(m,1,40000); q=reshape(n,1,40000); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% N_mcs=10; N=200; for mcs=1:N_mcs %%%How many Monte Carlo snnulahon (MCS) number parfor jj=1:N*N %%%%% N*N attempts are a MCS D2=D2{s(jj)}; D1=D1{s(jj)}; bb(jj)=D2(q(jj)); m1=s(jj); n1=q(jj); bb=bb(jj); if(n1==fix(n1/2)*2)==0 %due to two forms of hexagon. Determine the neighboring-matrix. NM1=[m1,n1-1;m1+1,n1;m1,n1+1;m1-1,n1+1;m1-1,n1;m1-1,n1-1]; else NM1=[m1+1,n1-1;m1+1,n1;m1+1,n1+1;m1,n1+1;m1-1,n1;m1,n1-1]; end NB2=[D2{NM1(1,1)} (NM1(1,2)), D2{NM1(2,1)}(NM1(2,2)), D2{NM1(3,1)}(NM1(3,2)),D2{NM1(4,1)}(NM1(4,2)), D2{NM1(5,1)}( NM1(5,2)), D2{NM1(6,1)}( NM1(6,2))];%% orientation number NB1=[D1{NM1(1,1)} (NM1(1,2)), D1{NM1(2,1)}(NM1(2,2)), D1{NM1(3,1)}(NM1(3,2)),D1{NM1(4,1)}(NM1(4,2)), D1{NM1(5,1)}( NM1(5,2)), D1{NM1(6,1)}( NM1(6,2))];%% different %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%