diff --git a/analysis.m b/analysis.m new file mode 100644 index 0000000000000000000000000000000000000000..e7eb1a2d40db6e60d44553f34eb1927067f86039 --- /dev/null +++ b/analysis.m @@ -0,0 +1,80 @@ + +%% distribution of cluster sizes N(n,t)= mean number of clusters of size n at time t +clear N; +N=zeros(700,size(c,2)); +T=1:size(c,2); +edges=0.5:1:700.5; +for I=T + for i=1:size(c,1) + N(:,I)=N(:,I)+histcounts(c{i,I},edges)'; + end +end + +I=find(sum(N,2),1,'last'); +N=N(1:I,:)/size(c,1); + + %% + figure(1) + clf; + I=2500; + n=1:size(N,1); + loglog(n,N(:,I)); + %image(N) + + + +%% mean number of clusters greater than k, number of free particles and size of largest cluster + +k=10; +figure(2) +clf; +subplot(3,1,1) + loglog(t,sum(N(k+1:end,:))) + hold on; + loglog(t(1500:2500),1000*(t(1500:2500)).^(-0.25)) + hold off; + +subplot(3,1,2) + semilogx(t,500-[2:size(N,1)]*N(2:end,:)) + ylim([0,500]) + +subplot(3,1,3) + + for i=T + largest(i)=find(N(:,i),1,'last'); + end + loglog(t,largest) + +%% mean cluster size bigger than k +k=10; +n=k+1:size(N,1); +tmp=sqrt(n/pi)*N(k+1:end,:); +%tmp=n*N(k+1:end,:); +R=tmp./sum(N(k+1:end,:)); + +figure(4) +loglog(t,R) +hold on +range=1600:length(t); +power=(4/15); +a=R(range(1))/t(range(1))^(power); +loglog(t(range),a*t(range).^(power)) +hold off +%xlim([1,1000]) +%ylim([1,10]) +xlabel('times (MCS)') +ylabel('linear cluster size') + +%% +%alternative method +%mean_size=cellfun(@(x) nanmean(x(x>k)),c); +%figure(5) +%loglog(t,nanmean(mean_size,1)) + +%xlim([1,100]) +%ylim([100,1000]) + +%% + + + \ No newline at end of file diff --git a/include/Particles.h b/include/Particles.h index 3b4a5a0a3e301a844e077aa603b787caf902904e..4e7ab2cbdb11a4e7eb22df329d753078c877cac5 100644 --- a/include/Particles.h +++ b/include/Particles.h @@ -129,6 +129,12 @@ public: return ori==1; } + + inline bool is_diffuse_ind(const int ind) const + { + return particles[ind]->ori==1; + } + inline bool is_bound(const int pos) const { int ori=0; @@ -139,8 +145,13 @@ public: ori = p-> ori; } return ori>1; + } + inline bool is_bound_ind(const int ind) const + { + return particles[ind]->ori>1; } + inline int get_orientation(const int pos) const { int ori=0; @@ -150,12 +161,24 @@ public: return ori-3; } + inline int get_orientation_ind(const int ind) const + { + return particles[ind]->ori-3; + } + inline void set_orientation(const int pos, const int ori) { assert(!grid1[pos].expired()); auto p(grid1[pos].lock()); p-> ori = ori; } + + inline void set_orientation_ind(const int ind, const int ori) + { + particles[ind]-> ori = ori; + } + + inline void set_pos(const int old_pos,const int new_pos, const int ind) { assert(particles[ind]->pos == old_pos); @@ -177,7 +200,7 @@ public: for(unsigned int i=0; i<n.size(); i++) { - if(int raw_ori=get_ori(n[i].position) && is_interaction_allowed(ori,raw_ori-3,n[i].slope)) + if(int raw_ori=get_ori(n[i].position); raw_ori>1 && is_interaction_allowed(ori,raw_ori-3,n[i].slope)) { count_bound++; } @@ -381,7 +404,7 @@ public: interactions.orientations.emplace_back(ori2); interactions.num_bonds++; interactions.num_diffuse++; - int bound_neigh_of_neigh = count_interacting_neighbors(n.position,ori); + int bound_neigh_of_neigh = count_interacting_neighbors(n.position,ori2); if(bound_neigh_of_neigh>1) { interactions.num_bonds=interactions.num_bonds+(bound_neigh_of_neigh-1); @@ -424,7 +447,7 @@ public: { // std::cout<<"binding"<<"\n"; binding_succ++; - set_orientation(pos,ori+3); + set_orientation_ind(ind,ori+3); for(unsigned int i=0; i<interactions.possible_interaction_pos.size(); i++) { set_orientation(interactions.possible_interaction_pos[i],interactions.orientations[i]+3); @@ -441,14 +464,13 @@ public: // get bound particle unbinding_attempt++; int particle_pos = get_pos(ind); - int ori = get_orientation(particle_pos); + int ori = get_orientation_ind(ind); Interactions interactions; interactions.num_bonds=0; interactions.num_diffuse=1; - int pos = get_pos(ind); - std::vector<Neighbour>& n(lattice.neighbours[pos]); + std::vector<Neighbour>& n(lattice.neighbours[particle_pos]); auto _is_bound = [this](const Neighbour &n) { @@ -479,7 +501,7 @@ public: if (rand<delta_E) { // std::cout<<"unbinding"<<"\n"; - set_orientation(particle_pos,1); + set_orientation_ind(ind,1); unbinding_succ++; if(interactions.possible_interaction_pos.size()>0) { @@ -492,9 +514,11 @@ public: } - void label(int ind, int i,std::vector<int> &labels,int &num_bonds) const + void label(int ind, int i,std::vector<int> &labels,int &num_bonds,int &num_particles) const { + labels[ind]=i; + num_particles++; int pos = get_pos(ind); int ori = get_orientation(pos); @@ -522,10 +546,11 @@ public: int inds = it-particles.begin(); return labels[inds]!=0; } + return false; }; - auto _label = [this,&labels,i,pos,&num_bonds](const Neighbour &n) + auto _label = [this,&labels,i,pos,&num_bonds,&num_particles](const Neighbour &n) { auto it=(std::find_if(begin(particles),end(particles), [n](std::shared_ptr<Particle> q) { @@ -537,19 +562,21 @@ public: int inds = it-particles.begin(); if(labels[inds]==0) { - label(inds,i,labels,num_bonds); + assert(particles[inds]->ori >1); + label(inds,i,labels,num_bonds,num_particles); } } }; - auto s1= n | std::views::filter(_is_bound) | std::views::filter(_is_allowed) | std::views::filter(_is_labelled); - auto cnt = std::ranges::distance(s1); + auto s1= n | std::views::filter(_is_bound) | std::views::filter(_is_allowed); + + auto s2=s1 | std::views::filter(_is_labelled); + auto cnt = std::ranges::distance(s2); num_bonds=num_bonds+cnt; - auto s2= n | std::views::filter(_is_bound) | std::views::filter(_is_allowed); - std::ranges::for_each(s2,_label); + std::ranges::for_each(s1,_label); } std::vector<int> orientations_vec() const diff --git a/main.cpp b/main.cpp index 473a53ecc763999030f68f7abbf33d5bb4831ede..aeae9a3c90ef0cbbea4faa1ae625dab724d62d25 100644 --- a/main.cpp +++ b/main.cpp @@ -31,14 +31,14 @@ int main(int argc,char *argv[]) } else { - alpha=0.1; J=2.6; + alpha=0.1; exponent=16; - slurm_index = 21; + slurm_index = 1001; std::cout << "Using default parameters." << '\n'; } density = 0.1; - const int MC_steps = 5*pow(10,4); // number of Monte Carlo Steps + const int MC_steps = 5*pow(10,5); // number of Monte Carlo Steps // const int MC_steps =500; int MC_counter = 0; // long double rand; @@ -63,13 +63,17 @@ int main(int argc,char *argv[]) // std::ofstream out2; // out2.open(fn2.str()); std::ostringstream fn; - fn << slurm_index << "_growth_dFrzB_bondspercluster_J_" << J << "_alpha_" << alpha << "_exp_" << exponent << ".txt";//k_un << "_" << k << ".txt"; + //fn << slurm_index << "_growth_dFrzB_bondspercluster_J_" << J << "_alpha_" << alpha << "_exp_" << exponent << ".txt";//k_un << "_" << k << ".txt"; + + fn << slurm_index << ".txt";//k_un << "_" << k << ".txt"; std::ofstream out; out.open(fn.str()); - std::ostringstream fn2; - fn2 << slurm_index << "_growth_dFrzB_labels_J_" << J << "_alpha_" << alpha << "_exp_" << exponent << ".txt";// k_un << "_" << k << ".txt"; - std::ofstream out2; - out2.open(fn2.str()); + +// std::ostringstream fn2; +// //fn2 << slurm_index << "_growth_dFrzB_labels_J_" << J << "_alpha_" << alpha << "_exp_" << exponent << ".txt";// k_un << "_" << k << ".txt"; +// fn2 << slurm_index << "_bonds.txt";//k_un << "_" << k << ".txt"; +// std::ofstream out2; +// out2.open(fn2.str()); // out << Nx << '\t' << Ny << '\t' <<'\n'; // out2 << Nx << '\t' << Ny << '\t' <<'\n'; @@ -127,14 +131,10 @@ int main(int argc,char *argv[]) // std::cout<<"ind \t"<<ind<<"\n"; rand=rand_size-ind; assert(rand<1); - if(particles.is_diffuse(particles.get_pos(ind))) + if(particles.is_diffuse_ind(ind)) { //Move diffusive particles - if(rand<0) - { - - } particles.attempt_diffusion(ind, rand); // print_container(particles.positions); @@ -142,16 +142,26 @@ int main(int argc,char *argv[]) //BINDING particles.attempt_binding(alpha,J,ind,rand); // print_container(particles.positions); + if(particles.is_bound(particles.get_pos(ind))){ - + assert(particles.get_orientation(particles.get_pos(ind))>-2); + assert(particles.get_ori(particles.get_pos(ind))==particles.particles[ind]->ori); + assert(particles.count_interacting_neighbors(particles.get_pos(ind),particles.get_orientation(particles.get_pos(ind)))); + } } - else if(particles.is_bound(particles.get_pos(ind))) + else if(particles.is_bound_ind(ind)) { // UNBINDING ATTEMPT + assert(particles.get_ori(particles.get_pos(ind))==particles.particles[ind]->ori); + assert(particles.count_interacting_neighbors(particles.get_pos(ind),particles.get_orientation(particles.get_pos(ind)))); particles.attempt_unbinding(alpha,J,ind,rand); // print_container(particles.positions); + if(particles.is_bound(particles.get_pos(ind))){ + assert(particles.get_ori(particles.get_pos(ind))==particles.particles[ind]->ori); + assert(particles.count_interacting_neighbors(particles.get_pos(ind),particles.get_orientation(particles.get_pos(ind)))); + } } @@ -168,31 +178,42 @@ int main(int argc,char *argv[]) - if(MC_counter%10000==0) + if(MC_counter%( MC_counter<100 ? 1: (int) pow(10,floor(log10(MC_counter))-2))==0) { - std::cout<<Nx<<"num of particles \t"<<particles.particles.size()<<'\n'; + + out << MC_counter << '\t'; + //out2 << MC_counter << '\t'; + + //std::cout<<Nx<<"num of particles \t"<<particles.particles.size()<<'\n'; // std::cout<<"Number of bonds: "; std::vector<int> labels(particles.particles.size(),0); int label_i = 1; for (unsigned int label_index=0; label_index < particles.particles.size(); label_index++) { - if(particles.is_bound(particles.get_pos(label_index)) && labels[label_index]==0) + if(particles.is_bound_ind(label_index) && labels[label_index]==0) { + assert(particles.is_bound(particles.get_pos(label_index))); + assert(particles.get_ori(particles.get_pos(label_index))==particles.particles[label_index]->ori); int num_bonds=0; - particles.label(label_index,label_i,labels,num_bonds); + int num_particles=0; + particles.label(label_index,label_i,labels,num_bonds,num_particles); label_i++; // std::cout << num_bonds << '\t'; - out << num_bonds<< '\t'; + assert(num_particles>1); + out << num_particles<< '\t'; + //out2 << num_bonds<< '\t'; } } //std::cout<<labels.size()<<" "<<particles.positions.size()<<'\n'; - std::cout<<"Number of clusters: " << std::ranges::max(labels) << '\n'; + //std::cout<<"Number of clusters: " << std::ranges::max(labels) << '\n'; // print_container(labels); out << '\n'; - particles.print_labels(out2,labels); - particles.print(out2); + //out2 << '\n'; + + //particles.print_labels(out2,labels); + //particles.print(out2); } MC_counter++; } diff --git a/main_script.m b/main_script.m new file mode 100644 index 0000000000000000000000000000000000000000..83516d3bb2b23b3d8960a173157529122932b897 --- /dev/null +++ b/main_script.m @@ -0,0 +1,37 @@ +clear; + + + + + +% parfor i=1:36 +% pause(0.0001*i) +% cmd=['set path=%path:C:\Program Files\MATLAB\R2020b\bin\win64;=% & bin\Release\frz_lattice_model.exe 2.6 0.1 16 ',num2str(i)]; +% system(cmd,'-echo');%outputs to current matlab directory +% end + + +%% +for i=1:1000 + j=1; + file=fopen([num2str(i),'.txt']); + tline=fgetl(file); + while ~isempty(tline) && ischar(tline) + z=textscan(tline,'%f','Delimiter','\t');%outputs cell array + t(j)=z{1}(1); + c{i,j}=z{1}(2:end); + tline=fgetl(file); + j=j+1; + end + +fclose(file); + +end + +save('results.mat','c','t'); + + +%% + + +