Untitled

#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <ctime>
#include <cstdlib>
using namespace std;

const int each_pictures = 100; //每幾步輸出一次 
const int x_unit = 100; //基板長 
const int y_unit = 100; //基板寬 
const int z_unit = 3; //基板高 
const int N = 10000; //模擬步數 
const int dep_layer = 2; //要進行磊晶的層(從0開始算)

int cubic[x_unit][y_unit][z_unit] = {}; //建立基板 
int pos[x_unit*y_unit*z_unit][3]={}; //存儲每個原子的位置 
int surface_pos[x_unit*y_unit][3]={}; //存儲磊晶層原子的位置
int event_list[x_unit*y_unit*z_unit*6][3][2]={}; //存儲所有可能發生的事件 
double k[x_unit*y_unit*z_unit*6]={}; //存儲所有發生事件的反應速率常數 
double p_list[x_unit*y_unit*z_unit*6]={}; //所有發生事件的機率列表 


int dep_times=0; //沉積過的次數 

string folder = "substrate/KMC_FCC/"; //資料夾位置 
string file="model_consequences_"; //檔名 

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

void build_A(int k)
{
	for(int i=0;i<x_unit;i++) 
		for(int j=0;j<y_unit;j++)
		{
			if(i%2==0)
			{
				if(j%6==0) cubic[i][j][k]=1;
			}
			else if(i%2==1)
			{
				if(j%6==3) cubic[i][j][k]=1;
			}
		}
}

void build_B(int k)
{
	for(int i=0;i<x_unit;i++) 
		for(int j=0;j<y_unit;j++)
		{
			if(i%2==0)
			{
				if(j%6==4) cubic[i][j][k]=1;
			}
			else if(i%2==1)
			{
				if(j%6==1) cubic[i][j][k]=1;
			}
		}
}

void build_C(int k)
{
	for(int i=0;i<x_unit;i++) 
		for(int j=0;j<y_unit;j++)
		{
			if(i%2==0)
			{
				if(j%6==2) cubic[i][j][k]=1;
			}
			else if(i%2==1)
			{
				if(j%6==5) cubic[i][j][k]=1;
			}
		}
}

void fill(int k)
{
	for(int i=0;i<x_unit;i++) 
		for(int j=0;j<y_unit;j++)
		{
			if(k<dep_layer&&cubic[i][j][k]==1) cubic[i][j][k]=2;
		}
}

void substrate_construction() //FCC基板建成 
{
	for(int k=0;k<z_unit;k++)
	{
		switch(k%3)
		{
			case 0:
				build_A(k);
				fill(k);
				build_B(k);
				break;
			case 1:
				build_B(k);
				fill(k);
				build_C(k);
				break;
			case 2:
				build_C(k);
				fill(k);
				build_A(k);
				break;
		}
	}
}

int counter() //回傳現有原子數 
{
	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) atoms++;
	return atoms;
}

int surface_counter() //回傳表面原子數 
{
	int atoms=0;
	for(int i=0;i<x_unit;i++)
		for(int j=0;j<y_unit;j++)
			if(cubic[i][j][dep_layer]==2) atoms++;
	return atoms;
}

void scan() //檢索所有原子並記錄在陣列 
{
	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)
				{
					pos[atoms][0]=i;
					pos[atoms][1]=j;
					pos[atoms][2]=k;
					atoms++;
				}
}

void surface_scan() //檢索表面原子並記錄在陣列 
{
	int atoms=0;
	for(int i=0;i<x_unit;i++)
		for(int j=0;j<y_unit;j++)
			if(cubic[i][j][dep_layer]==2)
			{
				surface_pos[atoms][0]=i;
				surface_pos[atoms][1]=j;
				surface_pos[atoms][2]=dep_layer;
				atoms++;
			}
}

double deposition_time() //定義沉積週期 
{
	double flux_rate = 10000;
	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(int storage) //建立機率列表 
{
	double probability[storage]={};
	double k_sum=0;
	
	for(int i=0;i<storage;i++)
		k_sum+=k[i];
		
	for(int i=0;i<storage;i++)
		probability[i]=k[i]/k_sum;
			
	for(int i=0;i<storage;i++)
	{
		p_list[i]=0; // 重置每次大迴圈的殘留資料 
		for(int j=0;j<i+1;j++)
		{
			p_list[i]+=probability[j];
		}
	}
	
	return k_sum;
}

int judgement(int storage) //判斷事件發生編號 
{
	float judge=float(rand()%1000)/float(1000);
	int tag=0; 
	for(int i=0;i<storage;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(int x,int y,int z)  //回傳位置上的原子狀態 0是不可存在原子 1是可存在原子 2是有原子 
{
	return cubic[x][y][z];
}

int bonding_numbers(int x,int y,int z,int n) //計算周圍鍵結數 [已考慮週期性] 
{
	int bn=0;
	
	int nn1[6][2] = {{ 1,-1}, {-2, 0}, { 1, 1}, {-1,-1}, { 2, 0}, {-1, 1}};
	int nn2[6][2] = {{-3,-1}, { 0,-2}, { 3,-1}, { 3, 1}, { 0, 2}, {-3, 1}};
	
	int c[6][2]={};
	
	if(n==1)
	{
		for(int i=0;i<6;i++)
		{
			c[i][0]=nn1[i][0];
			c[i][1]=nn1[i][1];
		}
	}
	else if(n==2)
	{
		for(int i=0;i<6;i++)
		{
			c[i][0]=nn2[i][0];
			c[i][1]=nn2[i][1];
		}
	}
	
	for(int i=0;i<6;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 event_nn1[][2], int event_nn2[][2], double Ea[][3]) //蒐集事件 
{
	int surface_atoms = surface_counter();
	int storage = 0;
		
	for(int i=0;i<surface_atoms;i++) //建立所有可能事件列表  
	{
		int ini_site[3] = {surface_pos[i][0],surface_pos[i][1],surface_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; //決定第N鄰近擴散
			
		for(int e=0;e<6;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]=surface_pos[i][0]+event_nn1[e][0];
				fin_site[1]=surface_pos[i][1]+event_nn1[e][1];
				periodic_site(fin_site);
			}
			else if(nn1_or_nn2==2) //第二鄰近擴散  
			{
				fin_site[0]=surface_pos[i][0]+event_nn2[e][0];
				fin_site[1]=surface_pos[i][1]+event_nn2[e][1];
				periodic_site(fin_site); 
			}
			
			int p_fin_site[3] = {fin_site[0],fin_site[1],surface_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(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, double t_dep, double k_sum, double t_loop) //事件執行 
{	
	if(t_loop>=t_dep*dep_times) //吸附事件發生 
	{
		dep_times++;
		
		int rx = rand()%x_unit;
		int ry = rand()%y_unit;
		
		while(cubic[rx][ry][dep_layer]!=1||bonding_numbers(rx,ry,dep_layer,1)>0) //骰子骰到磊晶層某個點還沒有原子為止 
		{
			rx = rand()%x_unit;
			ry = rand()%y_unit;
		}
		
		cubic[rx][ry][dep_layer]=2;
	}
	else //擴散事件發生 
	{
		t_loop+=1/k_sum;
		
		
		cubic[event_list[tag][0][0]][event_list[tag][1][0]][event_list[tag][2][0]]=0;	//移除原來的原子 
		cubic[event_list[tag][0][1]][event_list[tag][1][1]][event_list[tag][2][1]]=2;	//新增事件產生的原子 
	}
	
	return t_loop; //回傳經過時長(如果是沉積 經過時長不會變動) 
}

void output(int count, string folder_temp, string file_temp) //輸出檔案 
{
	int atoms = counter();
	scan();
	
	folder=folder_temp;
	file=file_temp;
	
	file+=trans(count);
	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 << "2 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 << "2 12.011" << endl << endl; // 1是基板原子 2是磊晶原子 
	outputFile << "Atoms # molecular" << endl << endl;
	
	for(int i=0;i<atoms;i++)
	{
		outputFile << i+1;
		if(pos[i][2]<dep_layer) outputFile << " 1 " << pos[i][0] << " " << pos[i][1] << " " << pos[i][2] << endl; // 磊晶層以上才會標註2 
		else outputFile << " 2 " << 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 event_nn1[6][2] = {{ 1,-1}, {-2, 0}, { 1, 1}, {-1,-1}, { 2, 0}, {-1, 1}}; //能移動的六個向量 
	int event_nn2[6][2] = {{-3,-1}, { 0,-2}, { 3,-1}, { 3, 1}, { 0, 2}, {-3, 1}}; //能移動的六個向量 

	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}};
                    
	
	int storage=0; //可能事件的數目 
	int tag=0; //發生事件的編號 
	double k_sum=0;
	double t_dep=0;
	double t_loop=0;
	
	t_dep = deposition_time(); //定義沉積週期 
	
	substrate_construction(); //FCC基板建成  
	
	for(int count=0;count<N;count++) //KMC循環開始 
	{
		surface_scan(); //紀錄表面原子位置 
			
		storage = event_collection(event_nn1,event_nn2,Ea); //蒐集所有可能發生的事件 
		
		k_sum = p_list_construction(storage); //建立所有事件的機率列表 
		
		tag = judgement(storage);  //決定發生哪一個擴散事件  
		
		t_loop = execution(r,tag,t_dep,k_sum,t_loop); //讓那個擴散事件發生 
		
		if(surface_counter()>((x_unit*y_unit*2)/5)) return 2; //第二停止條件 (當原子覆蓋率大於閾值時停止)
		
		if(count%each_pictures==0) output(count,folder_temp,file_temp); //每特定幾張做成.txt檔輸出 
	}
	
	return 1;