Untitled

mail@pastecode.io avatar
unknown
c_cpp
a month ago
7.7 kB
3
Indexable
Never
#include<bits/stdc++.h>
using namespace std;

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

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]==1) 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]==1)
				{
					p[atoms][0]=i;
					p[atoms][1]=j;
					p[atoms][2]=k;
					atoms++;
				}
}

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;

}

void 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];
		}
	}
}

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 bn=0;
	int c[2]={};
	
	c[0]=x+1, c[1]=y;
	periodic_site(c);
	if(cubic[c[0]][c[1]][z]==1) bn++;
	
	c[0]=x-1, c[1]=y;
	periodic_site(c);
	if(cubic[c[0]][c[1]][z]==1) bn++;
	
	c[0]=x, c[1]=y+1;
	periodic_site(c);
	if(cubic[c[0]][c[1]][z]==1) bn++;
	
	c[0]=x, c[1]=y-1;
	periodic_site(c);
	if(cubic[c[0]][c[1]][z]==1) 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[][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]}; //存取現有位置 
			
		for(int e=0;e<4;e++)
		{	 
			int fin_site[2]={pos[i][0]+event[e][0],pos[i][1]+event[e][1]};
			periodic_site(fin_site);
				
			int p_fin_site[3] = {fin_site[0],fin_site[1],pos[i][2]}; //存取可能位置 
				
			int ini_bn = bonding_numbers(ini_site[0],ini_site[1],ini_site[2]); //存取現有位置的鍵結數 
			int fin_bn = bonding_numbers(p_fin_site[0],p_fin_site[1],p_fin_site[2])-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])==1) 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; //回傳存儲事件數 
}

void execution(int r[], int tag, int event_list[][3][2]) //事件執行 
{
	float dep_or_dif = float(rand()%1000)/float(1000);
	
	if(dep_or_dif<0.02) //吸附事件發生 
	{
		int rx = rand()%x_unit;
		int ry = rand()%y_unit;
		
		while(cubic[rx][ry][0]==1) //骰子骰到那個點沒有原子為止 
		{
			rx = rand()%x_unit;
			ry = rand()%y_unit;
		}
		
		cubic[rx][ry][0]=1;
	}
	else //擴散事件發生 
	{
		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]]=1;
	}
}

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[4][2]={{1,1},{-1,-1},{-1,1},{1,-1}}; //能移動的四個向量 
	int event_list[x_unit*y_unit*z_unit*4][3][2]={}; //存儲所有可能發生的事件 
	int r_x = rand()%20+40;
	int r_y = rand()%20+40;
	int r[2] = {r_x,r_y};
	
	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; //發生事件的編號 
	
	cubic[r[0]][r[1]][0]=1;
	
	
	for(int count=0;count<N;count++) //KMC循環開始 
	{
		scan(pos); //紀錄所有原子位置 
			
		storage = event_collection(pos,event,Ea,k,event_list); //蒐集所有可能發生的事件 
		
		p_list_construction(k,p_list,storage); //建立所有事件的機率列表 
		
		tag = judgement(p_list,storage);  //決定發生哪一個事件  
		
		execution(r,tag,event_list); //讓那件事發生 
		
		if(count%10==0) output(pos,count,folder_temp,file_temp); //做成.txt檔輸出 
	}
	
	return 0;
}
Leave a Comment