Framework of my kMC¶
This is a framework of lipai’s kMC code
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<vector>
#include<time.h>
#define N_limit 80000000
using namespace std;
const int multiplier=50;
const int step_out=10000000;
const int row=2000;
const int col=2000;
const int c_num=10;
double t=0.0;
//const int G_row=200;
//const int G_col=200;
const double kT=0.112;
int num_evt[2]={0};
int mesh[row][col][2]={0};// save species info
int species_num[2]={0};
int num_diff[2]={0}; //save diffu events
double prob_each[3][4]={
//0 1 2
//diff C C2 e1
/* {0.00,0.00,0.00,0.00},
{0.20,1.00,0.00,0.00},//C
{0.00,0.00,0.00,2.30},//C2 e1:C2-> C+C
*/
{0.00, 0.00, 0.00, 0.00},
{0.50, 0.90, 0.00, 0.00},
{0.49, 0.00, 0.00, 2.80},
};
// C2->C+C 2.7eV, C+C->C2 ???
int sp_merge[2][2]={\
//1 2
//C C2
{2,0},//C
{0,0},//C2
};
int sp_e[2][2]={\
{0,0}, //C
{1,1}}; //C2
int evt_emerge[2][2]={\
//1 2
//C C2
{1,-1},//C
{-1,-1},//C2
};
int evt_spec[2]={-1,2};
int p[3][2];
int q[3][2];
void err(char* erro);
void init();
int get_nb(int x,int y,int xy,int i);
void find_evt();
void diffusion(int sp,int s_n,int evt);
void merge(int sp1,int s_n,int evt);
void spec_del(int sp,int s_n);
void creat(int sp,int x,int y);
void spec_e(int sp,int s_n,int evt);
void update_evt_nb(int x,int y);
void statistic(int n,int N);
int get_sp(int s_x,int s_y);
void add_c(int sp);
class species{
public:
species(int sp,int x,int y){
spec=sp;
s_x=x;
s_y=y;
prob[3]=prob_each[spec][3]; //desorp or decomp events
check_event();
}
void check_event(){
double r9;
for(int i=0;i<3;i++){ //3 directions
nb[i][0]=get_nb(s_x,s_y,0,i);
nb[i][1]=get_nb(s_x,s_y,1,i);
event[i]=get_sp(nb[i][0],nb[i][1]);
// if(event[i]==100) prob[i]=0; //2017 added for h
// else prob[i]=prob_each[spec][event[i]];
prob[i]=prob_each[spec][event[i]];
}
prob_s=0;
for(int i=0;i<4;i++) prob_s+=prob[i];
}
int update(){
int i;
double evt_sum=0;
double evt_add=0;
double evt_rand;
for(i=0;i<4;i++)
if(prob[i]!=0) evt_sum+=prob[i];
evt_rand=evt_sum*(double)rand()/RAND_MAX;
for(i=0;i<4;i++){
if(prob[i]!=0) evt_add+=prob[i];
if(evt_add>=evt_rand){
return i;
}
}
err("spec.update()");
return 0;
}
double prob[4],prob_s;
int event[3],nb[3][2];
int spec,s_x,s_y;
};
vector<vector<species> > spec_list(2);
int main(){
init();
int init_h=0,h;
for(int i=0;i<2;i++) init_h+=species_num[i];
for(int m=0;m<multiplier;m++){
for(int i=0;i<N_limit;i++){
// h=init_h;
find_evt();
// init_h=0;
// for(int j=0;j<2;j++) init_h+=species_num[j];
//if(init_h!=h) statistic(i,m);
if(i%step_out==0) statistic(i,0);
}
}
statistic(-1,0);
for(int i=0;i<2;i++){
printf("num_evt %d:\t%d\n",i,num_evt[i]);
printf("num_diff %d:\t%d\n",i,num_diff[i]);
}
return 0;
}
void err(char* erro){
printf("%s error!!!\n",erro);
scanf("%s",&erro);
}
void init(){
int i,j;
srand((unsigned)time(NULL));
for(i=1;i<3;i++) //2 species
for(j=0;j<4;j++)
if(prob_each[i][j]!=0)
prob_each[i][j]=2.71*pow(10,13)*exp(-prob_each[i][j]/kT);
for(i=1;i<3;i++)
for(j=1;j<3;j++)
if(prob_each[i][j]!=0)
prob_each[i][j]=prob_each[i][j]/2;
p[0][0]=-1;p[0][1]=1;// p for x even
p[1][0]=-1;p[1][1]=0;
p[2][0]=1;p[2][1]=0;
q[0][0]=-1;q[0][1]=0;//q for x odd
q[1][0]=1;q[1][1]=0;
q[2][0]=1;q[2][1]=-1;
for(i=0;i<c_num;i++) add_c(1);
//for(int i=row/2-G_row/2-1;i<row/2+G_row/2;i++)
//for(int j=col/2-G_col/2-1;j<col/2+G_col/2;j++)
//mesh[i][j][0]=9;
//mesh[row/2-G_row/2-1][col/2-G_col/2-1][0]=100;
//mesh[row/2-G_row/2-1][col/2+G_col/2][0]=100;
//mesh[row/2+G_row/2][col/2-G_col/2-1][0]=100;
//mesh[row/2+G_row/2][col/2+G_col/2][0]=100;
}
int get_nb(int x,int y,int xy,int i){
int x_temp,y_temp;
if(xy==0){
if(x%2==0) x_temp=x+p[i][0];
else x_temp=x+q[i][0];
if(x_temp>=row) return x_temp-=row;
else if(x_temp<0) return x_temp+=row;
else return x_temp;
}
else {
if(x%2==0) y_temp=y+p[i][1];
else y_temp=y+q[i][1];
if(y_temp>=col) return y_temp-=col;
else if(y_temp<0) return y_temp+=col;
else return y_temp;
}
}
void find_evt(){
double prob_rand;
double prob_sum=0;
double prob_add=0;
int drct,s_n;
double delta_t;
for(int i=0;i<2;i++)
if(!spec_list[i].empty())
for(int j=0;j<spec_list[i].size();j++) prob_sum+=spec_list[i][j].prob_s;
//printf("prob_s:\t%f\n",prob_sum);
//
do{
delta_t=(double)rand()/RAND_MAX;
}while(delta_t==0);
delta_t=-log(delta_t)/prob_sum;
t+=delta_t;
prob_rand=prob_sum*(double)rand()/RAND_MAX;
for(int i=0;i<2;i++)
if(!spec_list[i].empty()){
for(int j=0;j<spec_list[i].size();j++){
prob_add+=spec_list[i][j].prob_s;
if(prob_add>=prob_rand){
drct=spec_list[i][j].update();
if(drct<3){
if(spec_list[i][j].event[drct]==0){
diffusion(i+1,j,drct);
num_diff[i]++;
//printf("diffuse cxhy");
//
//printf("get it1\n");
}
else if(spec_list[i][j].event[drct]<=2){
merge(i+1,j,drct);
//printf("get it2\n");
//printf("merge cxhy\n");
}
else if(spec_list[i][j].event[drct]==100)
return;
else err("error in update1");
}
else {
spec_e(i+1,j,drct);
if(evt_spec[i]!=-1) num_evt[evt_spec[i]-1]++;
}
return;
}
}
}
//err("prob");
}
void diffusion(int sp,int s_n,int drct){
int x1,y1,x2,y2;
x1=spec_list[sp-1][s_n].s_x;
y1=spec_list[sp-1][s_n].s_y;
x2=spec_list[sp-1][s_n].nb[drct][0];
y2=spec_list[sp-1][s_n].nb[drct][1];
if(mesh[x1][y1][0]!=sp||mesh[x2][y2][0]!=0){ //for debug
//printf("in diffu,sp,x1,y1,x2,y2: %d\t%d\t%d\t%d\t%d\n",sp,x1,y1,x2,y2);
//if(mesh[x1][y1][0]!=sp) {printf("sp=%d,mesh=%d\n",sp,mesh[x1][y1][0]); err("diffusion1");}
if(mesh[x1][y1][0]!=sp) printf("dif1");
if(mesh[x2][y2][0]!=0) printf("dif1");
printf("drct=%d\tsp1=%d,mesh[0]=%d\n",drct,sp,mesh[x1][y1][0]);
printf("event[drct]:%d",spec_list[sp-1][s_n].event[drct]);
printf("sp2=%d,mesh2[1]=%d\n",mesh[x2][y2][0],mesh[x2][y2][1]);
printf("x2=%d,y2[1]=%d\n",spec_list[sp-1][mesh[x2][y2][1]].s_x,spec_list[sp-1][mesh[x2][y2][1]].s_y);
for(int i=0;i<3;i++){
printf("evt[i]=%d\n",spec_list[sp-1][s_n].event[i]);
update_evt_nb(x2,y2);
printf("evt[i]=%d\n",spec_list[sp-1][s_n].event[i]);
spec_list[sp-1][s_n].check_event();
printf("evt[i]=%d\n",spec_list[sp-1][s_n].event[i]);
printf("neigbour of 1: %d\t%d\n",get_nb(x1,y1,0,i),get_nb(x1,y1,1,i));
printf("neigbour of 2: %d\t%d\n",get_nb(x2,y2,0,i),get_nb(x2,y2,1,i));
}
err("diffusion");
//}
err("diffusion");
}
mesh[x2][y2][0]=mesh[x1][y1][0];
mesh[x2][y2][1]=mesh[x1][y1][1];
mesh[x1][y1][0]=0;
mesh[x1][y1][1]=0;
spec_list[sp-1][s_n].s_x=x2;
spec_list[sp-1][s_n].s_y=y2;
spec_list[sp-1][s_n].check_event();
update_evt_nb(x1,y1);
update_evt_nb(x2,y2);
}
void merge(int sp1,int s_n,int drct){
int sp2,sp_crt;
int x1,y1,x2,y2;
x1=spec_list[sp1-1][s_n].s_x;
y1=spec_list[sp1-1][s_n].s_y;
x2=spec_list[sp1-1][s_n].nb[drct][0];
y2=spec_list[sp1-1][s_n].nb[drct][1];
sp2=spec_list[sp1-1][s_n].event[drct];
spec_del(sp1,s_n);
spec_del(sp2,mesh[x2][y2][1]);
sp_crt=sp_merge[sp1-1][sp2-1];
if(evt_emerge[sp1-1][sp2-1]!=-1) num_evt[evt_emerge[sp1-1][sp2-1]-1]++;
if(sp_crt>0) creat(sp_crt,x1,y1);
else if(sp_crt==0) err("merge");
}
void spec_del(int sp,int s_n){
int x,y;
int spec_end;
species_num[sp-1]--;
x=spec_list[sp-1][s_n].s_x;
y=spec_list[sp-1][s_n].s_y;
if(sp!=mesh[x][y][0]) err("spec_del1");
if(s_n!=mesh[x][y][1]) err("spec_del2");
spec_end=spec_list[sp-1].size()-1;
if(s_n!=spec_end){
//mesh[spec_list[sp-1][spec_end].s_x][spec_list[sp-1][spec_end].s_y][0]=sp;
mesh[spec_list[sp-1][spec_end].s_x][spec_list[sp-1][spec_end].s_y][1]=s_n;
swap(spec_list[sp-1][s_n],spec_list[sp-1][spec_end]);
}
spec_list[sp-1].pop_back();
mesh[x][y][0]=0;
mesh[x][y][1]=0;
update_evt_nb(x,y);
}
void creat(int sp,int x,int y){
species_num[sp-1]++;
mesh[x][y][0]=sp;
mesh[x][y][1]=spec_list[sp-1].size();
spec_list[sp-1].push_back(species(sp,x,y));
update_evt_nb(x,y);
}
void spec_e(int sp,int s_n,int drct){
int sp1,sp2,xn,yn;
int x,y;
x=spec_list[sp-1][s_n].s_x;
y=spec_list[sp-1][s_n].s_y;
sp1=sp_e[sp-1][0];
sp2=sp_e[sp-1][1];
// revised 0925
if(sp1==0) err("spec_e");
spec_del(sp,s_n);
creat(sp1,x,y);
for(int i=0;i<3;i++){
xn=get_nb(x,y,0,i);
yn=get_nb(x,y,1,i);
if(mesh[xn][yn][0]==0){
creat(sp2,xn,yn);
break;
}
}
}
void update_evt_nb(int x,int y){
int sp,s_n;
int x_temp,y_temp;
for(int i=0;i<3;i++){
x_temp=get_nb(x,y,0,i);
y_temp=get_nb(x,y,1,i);
sp=mesh[x_temp][y_temp][0];
if(sp>0&&sp<=8){
s_n=mesh[x_temp][y_temp][1];
spec_list[sp-1][s_n].check_event();
}
}
}
void statistic(int n,int N){
int nn[2];
for(int i=0;i<2;i++) nn[i]=spec_list[i].size();
printf("time/c/c2:\t%e\t%d\t%d\n",t,species_num[0],species_num[1]);
}
int get_sp(int x,int y){
return mesh[x][y][0];
}
void add_c(int sp){
int x,y,xn,yn;
int flag=0;
do{
do{
x=(int)(row*(double)rand()/RAND_MAX);
y=(int)(col*(double)rand()/RAND_MAX);
}while(mesh[x][y][0]!=0);
for(int i=0;i<3;i++){
xn=get_nb(x,y,0,i);
yn=get_nb(x,y,1,i);
if(mesh[xn][yn][0]==0){
flag=1;
break;
}
}
}while(flag==0);
creat(sp,x,y);
}