Untitled

#include<bits/stdc++.h>
using namespace std;

const int x_unit = 50; //基板長 
const int y_unit = 50; //基板寬 
const int z_unit = 1; //基板高 
const int N = 100; //模擬次數 
string folder = "substrate/KMC_FCC/"; //資料夾位置 
string file="model_consequences_"; //檔名 
int cubic[x_unit][y_unit][z_unit] = {}; //建立基板 
int dep_times=0;

string trans(int num) //數字轉字串 
{
	string s = "";
	stringstream ss("");
	ss << num;
	ss >> s;
	return s;
}

int counter() //回傳現有原子數 
{
	int atoms=0;
	int p[x_unit*y_unit*z_unit][3]={};
	for(int i=0;i<x_unit;i++)
		for(int j=0;j<y_unit;j++)
			for(int k=0;k<z_unit;k++)
				if(cubic[i][j][k]==2) atoms++;
	return atoms;
}

void scan(int p[][3]) //檢索所有原子並記錄在陣列 
{
	int atoms=0;
	for(int i=0;i<x_unit;i++)
		for(int j=0;j<y_unit;j++)
			for(int k=0;k<z_unit;k++)
				if(cubic[i][j][k]==2)
				{
					p[atoms][0]=i;
					p[atoms][1]=j;
					p[atoms][2]=k;
					atoms++;
				}
}

double deposition_time()
{
	double flux_rate = 1000;
	double atoms_per_second = flux_rate*x_unit*y_unit;
	double t_dep = 1/atoms_per_second;
	
	return t_dep;
}

double Arrhenius_equation (double Ea, double T) //計算反應速率常數 
{
    double k = 0; //速率常數
    long long A = pow(10,12); //常數
    //double Ea = 0.5; //活化能 (ev)
    double R = 8.617*pow(10,-5); //氣體常數 (eV·K?1)
    //double T = 298; //絕對溫度 (k)

    k = A * exp((-Ea)/(R*T));

    return k;

}

double p_list_construction(double k[],double p_list[],int atoms) //建立機率列表 
{
	double probability[atoms*4]={};
	double k_sum=0;
	
	for(int i=0;i<atoms*4;i++)
		k_sum+=k[i];
		
	for(int i=0;i<atoms*4;i++)
		probability[i]=k[i]/k_sum;
			
	for(int i=0;i<atoms*4;i++)
	{
		p_list[i]=0; // 重置每次大迴圈的殘留資料 
		for(int j=0;j<i+1;j++)
		{
			p_list[i]+=probability[j];
		}
	}
	
	return k_sum;
}

int judgement(double p_list[],int atoms) //判斷事件發生編號 
{
	float judge=float(rand()%1000)/float(1000);
	int tag=0; 
	for(int i=0;i<atoms*4;i++)
		{
			tag=i;
			if(judge<p_list[i]) break;
		}
	return tag;
}

void periodic_site(int r[]) //產生週期性(需放入一維[2]陣列) 
{
	r[0]%=x_unit;
	r[1]%=y_unit;
	if(r[0]<0) r[0]+=x_unit;
	if(r[1]<0) r[1]+=y_unit;
}

int this_site_is_occupied(int x,int y,int z)  //回傳位置上是否有原子 
{
	return cubic[x][y][z];
}

int bonding_numbers(int x,int y,int z,int n) //計算周圍(前、後、左、右)鍵結數 [已考慮週期性] 
{
	int bn=0;
	
	int c[4][2]={{n,0},{-n,0},{0,n},{0,-n}};
	
	for(int i=0;i<4;i++)
	{
		c[i][0]+=x;
		c[i][1]+=y;
		
		int p_c[2]={c[i][0],c[i][1]};
		periodic_site(p_c);
		
		if(cubic[p_c[0]][p_c[1]][z]==2) bn++;
	}
	
	return bn;
}

double Ea_we_need(double Ea[16][3],int initial_bn,int final_bn) //選取對應能障 
{
	for(int i=0;i<16;i++)
		if(Ea[i][0]==initial_bn&&Ea[i][1]==final_bn) return Ea[i][2];
		
}

int event_collection(int pos[][3], int event_nn1[][2], int event_nn2[][2], double Ea[][3], double k[], int event_list[][3][2]) //蒐集事件 
{
	int atoms=counter();
	int storage = 0;
		
	for(int i=0;i<atoms;i++) //建立所有可能事件列表  
	{
		int ini_site[3] = {pos[i][0],pos[i][1],pos[i][2]}; //存取現有位置 
		int nn1_or_nn2=1;
		if(bonding_numbers(ini_site[0],ini_site[1],ini_site[2],2)>0) nn1_or_nn2=2;
			
		for(int e=0;e<4;e++)
		{	 
			int ini_bn = bonding_numbers(ini_site[0],ini_site[1],ini_site[2],nn1_or_nn2); //存取現有位置的鍵結數
			
			int fin_site[2]={};
			
			if(nn1_or_nn2==1)
			{
				fin_site[0]=pos[i][0]+event_nn1[e][0];
				fin_site[1]=pos[i][1]+event_nn1[e][1];
				periodic_site(fin_site);
			}
			else if(nn1_or_nn2==2)
			{
				fin_site[0]=pos[i][0]+event_nn2[e][0];
				fin_site[1]=pos[i][1]+event_nn2[e][1];
				periodic_site(fin_site); 
			}
			
			int p_fin_site[3] = {fin_site[0],fin_site[1],pos[i][2]}; //存取可能位置 
				
			int fin_bn = bonding_numbers(p_fin_site[0],p_fin_site[1],p_fin_site[2],nn1_or_nn2)-1; //存取可能位置的鍵結數
				
			int k_temp = Arrhenius_equation(Ea_we_need(Ea,ini_bn,fin_bn),298);
			
			if(this_site_is_occupied(p_fin_site[0],p_fin_site[1],p_fin_site[2])==2) k_temp=0;
			if(ini_bn>3) k_temp=0;
				
			if(k_temp!=0) //存取所有可能事件及其反應速率常數  //event_list[][3][2]: [x][y][z] ([][][0]), [pb_x][pb_y][pb_z] ([][][1]) 
			{
				k[storage]=k_temp;
				event_list[storage][0][0]= ini_site[0];
				event_list[storage][1][0]= ini_site[1];
				event_list[storage][2][0]= ini_site[2];
				event_list[storage][0][1]= p_fin_site[0];
				event_list[storage][1][1]= p_fin_site[1];
				event_list[storage][2][1]= p_fin_site[2];
				storage++;
			}
		}
	}
	
	return storage; //回傳存儲事件數 
}

double execution(int r[], int tag, int event_list[][3][2], double t_dep, double k_sum, double t_loop) //事件執行 
{	
	if(t_loop>=t_dep*dep_times) //吸附事件發生 
	{
		dep_times++;
		
		cout << "dep ";
		
		int rx = rand()%x_unit;
		int ry = rand()%y_unit;
		
		while(cubic[rx][ry][0]==2) //骰子骰到那個點沒有原子為止 
		{
			rx = rand()%x_unit;
			ry = rand()%y_unit;
		}
		
		cubic[rx][ry][0]=2;
	}
	else //擴散事件發生 
	{
		t_loop+=1/k_sum;
		
		cout << "dif " << 1/k_sum << " ";
		
		cubic[event_list[tag][0][0]][event_list[tag][1][0]][event_list[tag][2][0]]=0;
		
		r[0]=event_list[tag][0][1];
    	r[1]=event_list[tag][1][1];
    	periodic_site(r);
    	
    	cubic[r[0]][r[1]][event_list[tag][2][1]]=2;
	}
	
	return t_loop;
}

void output(int pos[][3], int c, string folder_temp, string file_temp) //輸出檔案 
{
	int atoms=counter();
	scan(pos);
	
	folder=folder_temp;
	file=file_temp;
	
	file+=trans(c);
	file+=".txt";
	folder+=file;
	
	ofstream outputFile(folder.c_str());
	outputFile << "LAMMPS data file via write_data, version 23 Jun 2022, timestep = 0" << endl << endl;
	
	outputFile << atoms << " atoms" << endl << "1 atom types" << endl << endl;
	
	outputFile << "0 " << x_unit-1 << " xlo xhi" << endl;
	outputFile << "0 " << y_unit-1 << " ylo yhi" << endl;
	outputFile << "0 " << z_unit-1 << " zlo zhi" << endl << endl;
	outputFile << "Masses" << endl << endl << "1 12.011" << endl << endl;
	outputFile << "Atoms # molecular" << endl << endl;
	
	for(int i=0;i<atoms;i++)
		outputFile << i+1 << " 1 " << pos[i][0] << " " << pos[i][1] << " " << pos[i][2] << endl;
	
	outputFile.close();
}

int main()
{
	srand(time(NULL));
	
	string folder_temp = folder;//備份資料夾名稱 
	string file_temp = file;//備份檔案名稱 
	int pos[x_unit*y_unit*z_unit][3]={}; //存儲每個原子的位置 
	int event_nn1[4][2]={{1,0},{-1,0},{0,1},{0,-1}}; //能移動的四個向量 
	int event_nn2[4][2]={{2,0},{-2,0},{0,2},{0,-2}}; //能移動的四個向量
	int event_list[x_unit*y_unit*z_unit*4][3][2]={}; //存儲所有可能發生的事件 
//	int r_x = rand()%((x_unit*1)/5)+((x_unit*2)/5);
//	int r_y = rand()%((y_unit*1)/5)+((y_unit*2)/5);
	int r[2] = {};
	
	double Ea[16][3]={{0,0,0.29}, {0,1,0.29}, {0,2,0.29}, {0,3,0.29}, //不同鍵結數下的能障大小  // {初始鍵結數,結束鍵結數,能障}
                      {1,0,0.40}, {1,1,0.35}, {1,2,0.33}, {1,3,0.30},
                      {2,0,0.50}, {2,1,0.45}, {2,2,0.43}, {2,3,0.40},
                      {3,0,0.60}, {3,1,0.58}, {3,2,0.55}, {3,3,0.53}};
                    
	double k[x_unit*y_unit*z_unit*4]={}; //存儲所有發生事件的反應速率常數 
	double p_list[x_unit*y_unit*z_unit*4]={}; //所有發生事件的機率列表 
	int storage=0; //可能事件的數目 
	int tag=0; //發生事件的編號 
	double k_sum=0;
	double t_dep=0;
	int dep_times=1;
	double t_loop=0;
	
	//cubic[r[0]][r[1]][0]=2;
	
	t_dep = deposition_time();
	
	
	for(int count=0;count<N;count++) //KMC循環開始 
	{
		scan(pos); //紀錄所有原子位置 
			
		storage = event_collection(pos,event_nn1,event_nn2,Ea,k,event_list); //蒐集所有可能發生的事件 
		
		k_sum = p_list_construction(k,p_list,storage); //建立所有事件的機率列表 
		
		tag = judgement(p_list,storage);  //決定發生哪一個事件  
		
		t_loop = execution(r,tag,event_list,t_dep,k_sum,t_loop); //讓那件事發生 
		
		cout << "After "<< count << " loop:" << " // " << storage << " events" << " // " << " t_loop: " << t_loop << endl;
		
		output(pos,count,folder_temp,file_temp); //做成.txt檔輸出 
	}
	return 0;
}