Untitled
unknown
c_cpp
20 days ago
8.7 kB
4
Indexable
Never
#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; }
Leave a Comment