Framework of my kMC
===========

This is a framework of my 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);
      }