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