⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 111.txt

📁 晶体生长蒙特卡洛模拟机制 希望和大家共享
💻 TXT
📖 第 1 页 / 共 3 页
字号:
🔍
123 
	int	a,b,c,I,J;
	unsigned char atom;

	for(i=0;i<sizex;i++) for(j=0;j<sizey;j++){
		atom=config(&a,&b,&c,i,j);//如果底部是A,c=0;B，如果是，16
		k=a*(11-a)/2+b;
		if(atom=='A'){
			k+=c;
		} else {
			k+=32+c;
		}
		Config[i][j]=k;
	}
}
surface_init()//SPI
{

	int	i,j,k,sz;
unsigned char	A,B,T;
unsigned char	*pt;

	for(i=BASE;i<HORIZON;i++) for(j=0;j<sizex;j++) for(k=0;k<sizey;k++)
	{
		Crystal[i][j][k]=0;
	}
	A='A';B='B';
        i=BASE;
        for(j=0;j<sizex;j++){for(k=0;k<sizey;k++){
			if( INIT1 > R_RANDOM()*100.0 ) Crystal[i][j][k]=B;
                                else    Crystal[i][j][k]=A;
        }}
        for(i=0;i<=BASE-1;i++){for(j=0;j<sizex;j++){for(k=0;k<sizey;k++){
			if( INIT1 > R_RANDOM()*100.0 ) Crystal[i][j][k]=B;
                                else    Crystal[i][j][k]=A;
        }}}
        for(j=0;j<sizex;j++) for(k=0;k<sizey;k++) {
                Surface[j][k]=BASE;
	}
}
//Surface[x][y]*(4=Bius)+Base_h[x][y] 如此它测试高度。
/*周期性边界条件和集合为介绍步*/
Spiral_Base_Set(){
	int m,n,h,hh,i,k;
	m=sizex/2;n=sizey/2;
	h=0;
	for(i=Step_N;i>0;i--){	//Spiral source
		hh=h;
		for(k=0;k<i;k++){Base_h[m-k][n-i+1]=hh;hh++;}
		for(k=0;k<i;k++){Base_h[m-i+1][n-i+2+k]=hh;hh++;}
		for(k=0;k<i;k++){Base_h[m-i+2+k][n+1]=hh;hh++;}
		for(k=0;k<i;k++){Base_h[m+1][n-k]=hh;hh++;}
		h=h+3;
	}
	for(i=0;i<n-Step_N+1;i++)Base_h[m+1][i]=Step_N*4-1;
	for(i=0;i<=n+1;i++){			//右上整形
		k=m+2;
		while(Base_h[k-1][i]>0){Base_h[k][i]=Base_h[k-1][i]-1;k++;}
	}
	for(i=m-1;i>=m-Step_N+1;i--){	//上中間整形
		k=n-Step_N;
		while(Base_h[i][k-1]>0){Base_h[i][k]=Base_h[i][k-1]-1;k--;}
	}
	for(i=0;i<=n+1;i++){
		k=m-Step_N;
		while(Base_h[k+1][i]>0){Base_h[k][i]=Base_h[k+1][i]-1;k--;}
	}
	for(i=0;i<sizex;i++){
		k=n+2;
		while(Base_h[i][k-1]>0){Base_h[i][k]=Base_h[i][k-1]-1;k++;}
	}

}
Step_Base_Set(){
	int i,j,k,m,jj,n;
	k=sizex/Step_N;n=Step_N;
	jj=0;
	for(i=k/2;i<sizex;i += k){
		for(j=jj;j<=i;j++)Color_Base[j]=n;
		n--;jj=i+1;
		for(m=1;m<=3;m++){
			for(j=0;j<sizey;j++){
				Base_h[i-3+m][j]=m;
			}
		}
	}
	for(j=i;j<sizex;j++)Color_Base[j]=n;
}
perio_set()
{
	int	i,k;
	
	for(i=1;i<=sizex;i++) P_X[i]=i-1;for(i=1;i<=sizey;i++)P_Y[i]=i-1;
	for(i=0;i<sizex;i++)for(k=0;k<sizey;k++){Base_h[i][k]=0;}
	for(i=0;i<sizex;i++) Color_Base[i]=0;
	if(MODEL==1){
		SPI_UP_x=sizex-2;SPI_LOW_x=1;	SPI_UP_y=sizey-2;SPI_LOW_y=1;
		P_X[0]=1;P_X[sizex+1]=sizex-2;	P_Y[0]=1;P_Y[sizey+1]=sizey-2;//Spiral; 对称界限
		Spiral_Base_Set();
	} else {
		SPI_UP_x=sizex;SPI_LOW_x=-1;	SPI_UP_y=sizey;SPI_LOW_y=-1;
		P_X[0]=sizex-1;P_X[sizex+1]=0;	P_Y[0]=sizey-1;P_Y[sizey+1]=0;//周期境界
		if(Step_N>0){
			Step_Base_Set();
		}
	}
}
int furoku;
prob_def()
{
	int	i,j,k;
	double ramda_A,ramda_B,kmxa,kmxb;
	
	k_a=(1-alpha)*F;k_b=alpha*F;
	ramda_A=ramda_a*ramda_a;
	ramda_B=ramda_b*ramda_b;
	for(k=0;k<STATE;k++) {
		Kxij[k]=0.0;choice_flg[k]=Prob[k][0]=Prob[k][1]=Prob[k][2]=0;
	}
	for(i=0;i<5;i++)for(j=0;j<5;j++){
		if(i+j>4) continue;		//这里4。1它下
		k=i*(11-i)/2+j;		//这个系统是重要的
		choice_flg[k]=choice_flg[k+16]=choice_flg[k+32]=choice_flg[k+48]=1;
		kmxa=exp(-i*fai_aa-j*fai_ab);kmxb=exp(-i*fai_ab-j*fai_bb);
	Kxij[k]=k_a+k_b+(1.0+ramda_A)*kmxa*exp(-fai_aa);// 在A A之下
	Kxij[k+16]=k_a+k_b+(1.0+ramda_A)*kmxa*exp(-fai_ab);//在A B之下
		Kxij[k+32]=k_a+k_b+(1.0+ramda_B)*kmxb*exp(-fai_ab_v);// 在A之下Kxij[k+48]=k_a+k_b+(1.0+ramda_B)*kmxb*exp(-fai_bb);// 在B之下
	Prob[k][0]=ramda_A*kmxa*exp(-fai_aa)/Kxij[k];// 表面传播	
	Prob[k+16][0]=ramda_A*kmxa*exp(-fai_ab)/Kxij[k+16];
		Prob[k+32][0]=ramda_B*kmxb*exp(-fai_ab_v)/Kxij[k+32];
		Prob[k+48][0]=ramda_B*kmxb*exp(-fai_bb)/Kxij[k+48];
	Prob[k][1]=Prob[k][0]+k_a/Kxij[k];// 可能性A原子落的地方
	Prob[k][2]=Prob[k][1]+k_b/Kxij[k];//可能性B原子落的地方，除那蒸发之外
		Prob[k+16][1]=Prob[k+16][0]+k_a/Kxij[k+16];
		Prob[k+16][2]=Prob[k+16][1]+k_b/Kxij[k+16];
		Prob[k+32][1]=Prob[k+32][0]+k_a/Kxij[k+32];
		Prob[k+32][2]=Prob[k+32][1]+k_b/Kxij[k+32];
		Prob[k+48][1]=Prob[k+48][0]+k_a/Kxij[k+48];
		Prob[k+48][2]=Prob[k+48][1]+k_b/Kxij[k+48];
	}
}
static int c_ount=0;
/*模仿物质 */
monte()
{
	int		k,member,I,J;
	double r;

for(c_ount=0;c_ount<Limit;c_ount++){
	Time += 1.0/State_tree[1];
	k=choice();
	member=M_RANDOM(Nxij[k])+Set_pt[k];
	I=Surface_tblX[member];	J=Surface_tblY[member];
	r=R_RANDOM();
//printf("%e %e k=%d\n",r,Prob[k][2],k);getch();
	if( r < Prob[k][0] ){
		diffusion(I,J,k);
	} else {
		if( r <= Prob[k][1] ){
			create_A(I,J); 	C_a++;
		} else {
			if( r <= Prob[k][2] ) {
				create_B(I,J);C_b++;
			} else {
				if( Crystal[Surface[I][J]][I][J]=='A' ) {
					E_a++;
				} else {
					E_b++;
				}
				evapolation(I,J);
			}
		}
	}//getch();
	if(timing==1){
		for(k=1;k<Speed;k++)INKEY();
	}
}
}
static	int	x[]={0,0,1,-1,1,-1,1,-1};
static	int	y[]={1,-1,0,0,1,1,-1,-1};
evapolation(I,J)
	int	I,J;
{
	int	a,b,c,i,j;
	unsigned char	atom,old_config,new_config;
	int		new_hight,hight,hbius;
	int	(*action)(int,int,int),Del_A(),Del_B();

	new_hight=--Surface[I][J];
        if( new_hight == 0 ) { printf("Under flow Error\n");EXIT(1);}
	hight=new_hight+1;
	hbius=hight*Bius;
	atom=Crystal[hight][I][J];
	Crystal[hight][I][J] = 0;
	if( atom == 'A' ) action=Del_A; else action=Del_B;
	
	i=P_X[I+2];
	if( ((Surface[i][J]*Bius+Base_h[i][J])-(hbius+Base_h[I][J]))/2==0 ){
		if(I<SPI_UP_x){action(i,J,0);} else {if(I==SPI_UP_x)action(i,J,1);}
	}

	i=P_X[I];
	if( ((Surface[i][J]*Bius+Base_h[i][J])-(hbius+Base_h[I][J]))/2==0 ){
		if(I>SPI_LOW_x){action(i,J,0);} else {if(I==SPI_LOW_x)action(i,J,1);}
	}

	j=P_Y[J+2];
	if( ((Surface[I][j]*Bius+Base_h[I][j])-(hbius+Base_h[I][J]))/2==0 ){
		if(J<SPI_UP_y){action(I,j,0);} else {if(J==SPI_UP_y)action(I,j,1);}
	}

	j=P_Y[J];
	if( ((Surface[I][j]*Bius+Base_h[I][j])-(hbius+Base_h[I][J]))/2==0 ){
		if(J>SPI_LOW_y){action(I,j,0);} else {if(J==SPI_LOW_y)action(I,j,1);}
	}

	old_config=Config[I][J];
		i=G[old_config&15][0]+G[old_config&15][1];
#ifdef	DEBUG
if(i<0 || i>4)printf("??????????????evapo %d\n",i);
#endif
		kink_E[i]++;
	atom=config(&a,&b,&c,I,J);
	if( atom == 'A' ) {
		new_config=a*(11-a)/2+b+c;
		atom=A_atom;
	} else {
		new_config=a*(11-a)/2+b+32+c;
		atom=B_atom;
	}
#ifdef	DEBUG
if(choice_flg[new_config]==0)printf("Error in eva,%d\n",new_config);
#endif
#ifdef	GR
	if(timing==1){
		PSET_BE(I,J,new_hight+Color_Base[I],atom );
        PSETW(I,J,(new_hight+CNP+Color_Base[I])%CN+atom );
	}
#endif
	if( old_config != new_config ) modify(old_config,new_config,I,J);
}
Del_A(I,J,sw)
	int	I,J,sw;
{
	int	i,j;
	unsigned char old_config,new_config;
	
	old_config = Config[I][J];
	i=G[old_config&15][0];j=G[old_config&15][1];
	new_config=old_config-i*(11-i)/2+(i-1-sw)*(12-i+sw)/2;
#ifdef	DEBUG
if(choice_flg[new_config]==0)printf("Error in dela,%d\n",new_config);
#endif
	modify(old_config,new_config,I,J);
}
Del_B(I,J,sw)
	int	I,J,sw;
{
	int	i,j;
	unsigned char old_config,new_config;
	
	old_config = Config[I][J];
	new_config=old_config-1-sw;
#ifdef	DEBUG
if(choice_flg[new_config]==0)printf("Error in delb,%d\n",new_config);
#endif
	modify(old_config,new_config,I,J);
}
Reset(){int i,j,k,s_min,s_max;
	s_min=s_max=Surface[0][0];
	for(i=0;i<sizex;i++)for(j=0;j<sizey;j++){
		if(s_min>Surface[i][j])s_min=Surface[i][j];
		if(s_max<Surface[i][j])s_max=Surface[i][j];
	}
	if(s_min<=BASE){printf("Can not reset the surface\n");EXIT(2);}
	for(k=0;k<=s_max-s_min+BASE;k++)
	for(i=0;i<sizex;i++)for(j=0;j<sizey;j++){
		Crystal[k][i][j]=Crystal[k+s_min-BASE][i][j];
	}
	for(i=0;i<sizex;i++)for(j=0;j<sizey;j++){
		Surface[i][j] -=s_min-BASE;
	}
	for(k=s_max-s_min+BASE+1;k<HORIZON;k++)
	for(i=0;i<sizex;i++)for(j=0;j<sizey;j++){
		Crystal[k][i][j]=0;
	}
}
create_A(I,J)
	int	I,J;
{
	int	a,b,c,i,j,new_hight,hbius;
	unsigned char	atom,old_config,new_config;
	
	if( Surface[I][J] == HORIZON-1 ) {Reset();}
	new_hight=++Surface[I][J];
	hbius=new_hight*Bius;
	Crystal[new_hight][I][J] = 'A';

	i=P_X[I+2];
	if( ((Surface[i][J]*Bius+Base_h[i][J])-(hbius+Base_h[I][J]))/2==0 ){
		if(I<SPI_UP_x){Add_A(i,J,0);} else {if(I==SPI_UP_x)Add_A(i,J,1);}
	}

	i=P_X[I];
	if( ((Surface[i][J]*Bius+Base_h[i][J])-(hbius+Base_h[I][J]))/2==0 ){
		if(I>SPI_LOW_x){Add_A(i,J,0);} else {if(I==SPI_LOW_x)Add_A(i,J,1);}
	}

	j=P_Y[J+2];
	if(((Surface[I][j]*Bius+Base_h[I][j])-(hbius+Base_h[I][J]))/2==0){
		if(J<SPI_UP_y){Add_A(I,j,0);} else {if(J==SPI_UP_y) Add_A(I,j,1);}
	}

	j=P_Y[J];
	if(((Surface[I][j]*Bius+Base_h[I][j])-(hbius+Base_h[I][J]))/2==0){
		if(J>SPI_LOW_y){Add_A(I,j,0);} else {if(J==SPI_LOW_y)Add_A(I,j,1);}
	}
	old_config=Config[I][J];
	atom=config(&a,&b,&c,I,J);
	new_config=a*(11-a)/2+b+c;
	i=G[new_config&15][0]+G[new_config&15][1];
#ifdef	DEBUG
		if(i<0 || i>4)printf("??????????????C_A %d\n",i);
		if(choice_flg[new_config]==0)printf("Error in crea,%d\n",new_config);
#endif
	kink_C[i]++;
#ifdef	GR
	if(timing==1){
		PSET_BH(I,J,new_hight+Color_Base[I],A_atom);
		PSETW(I,J,(new_hight+CNP+Color_Base[I])%CN+A_atom);
	}
#endif
	if( old_config != new_config ) modify(old_config,new_config,I,J);
}
create_B(I,J)
	int	I,J;
{
	int	a,b,c,i,j,new_hight,hbius;
	unsigned char	atom,old_config,new_config;
	if( Surface[I][J] == HORIZON-1 ) {Reset();}
	new_hight=++Surface[I][J];
	hbius=new_hight*Bius;
	Crystal[new_hight][I][J] = 'B';

	i=P_X[I+2];
	if( ((Surface[i][J]*Bius+Base_h[i][J])-(hbius+Base_h[I][J]))/2==0 ){
		if(I<SPI_UP_x){Add_B(i,J,0);} else {if(I==SPI_UP_x)Add_B(i,J,1);}
	}

	i=P_X[I];
	if( ((Surface[i][J]*Bius+Base_h[i][J])-(hbius+Base_h[I][J]))/2==0 ){
		if(I>SPI_LOW_x){Add_B(i,J,0);} else {if(I==SPI_LOW_x)Add_B(i,J,1);}
	}

	j=P_Y[J+2];
	if(((Surface[I][j]*Bius+Base_h[I][j])-(hbius+Base_h[I][J]))/2==0){
		if(J<SPI_UP_y){ Add_B(I,j,0);} else {if(J==SPI_UP_y)Add_B(I,j,1);}
	}

	j=P_Y[J];
	if(((Surface[I][j]*Bius+Base_h[I][j])-(hbius+Base_h[I][J]))/2==0){
		if(J>SPI_LOW_y){ Add_B(I,j,0);} else {if(J==SPI_LOW_y)Add_B(I,j,1);}
	}
	old_config=Config[I][J];
	atom=config(&a,&b,&c,I,J);
	new_config=a*(11-a)/2+b+c+32;
	i=G[new_config&15][0]+G[new_config&15][1];
#ifdef	DEBUG
		if(i<0 || i>4)printf("?????????????? C_B %d\n",i);
		if(choice_flg[new_config]==0)printf("Error in creb,%d\n",new_config);
#endif

	kink_C[i]++;

#ifdef	GR
	if(timing==1){
		PSET_BH(I,J,new_hight+Color_Base[I],B_atom);
		PSETW(I,J,(new_hight+CNP+Color_Base[I])%CN+B_atom);
	}
#endif

	if( old_config != new_config ) modify(old_config,new_config,I,J);
}
Add_A(I,J,sw)
	int	I,J,sw;
{
	int	i,j;
	unsigned char old_config,new_config;
	
	old_config = Config[I][J];
	i=G[old_config&15][0], j=G[old_config&15][1];
	new_config=old_config-i*(11-i)/2+(i+1+sw)*(10-i-sw)/2;

#ifdef	DEBUG
if(choice_flg[new_config]==0)printf("Error in adda,%d i=%d j=%d I=%d J=%d sw=%d\n",
											new_config,i,j,I,J,sw);
#endif
	modify(old_config,new_config,I,J);
}
Add_B(I,J,sw)
	int	I,J,sw;
{
	int	i,j;
	unsigned char old_config,new_config;
	
	old_config = Config[I][J];
	new_config=old_config+1+sw;
#ifdef	DEBUG
if(choice_flg[new_config]==0)printf("Error in addB,%d I=%d J=%d sw=%d\n",
											new_config,I,J,sw);
#endif
	modify(old_config,new_config,I,J);
}
modify(old_config,new_config,I,J)
	unsigned char	old_config,new_config;
	int	I,J;
{
	double	hindo;
	int	k,st,n,i,j;

	
	k=old_config;
	hindo=Kxij[k];
	k+=64;
	while( k>0 ){
		State_tree[k] -= hindo;
		k /= 2;
	}
	k=new_config;
	hindo=Kxij[k];
	k+=64;
#ifdef	DEBUG
	if( hindo == 0.0 ){
		printf("Fatal Error in modify\n");
		printf("New=%d Old=%d I=%d J=%d\n",new_config,old_config,I,J);
		printf("O %e N %e\n",Kxij[old_config],Kxij[new_config]);
		exit(1);
	}
#endif
	while( k>0 ){
		State_tree[k] += hindo;
		k /= 2;
	}
	Config[I][J] = new_config;
	n=--Nxij[old_config]+Set_pt[old_config];	/*删除老 */
	st=Surface_pt[I][J];
	Surface_tblX[st]=i=Surface_tblX[n];
	Surface_tblY[st]=j=Surface_tblY[n];
	Surface_pt[i][j]=st;

	n=Nxij[new_config]++;				/*增加新 */
	st=Set_pt[new_config]+n;
	if( st>=Set_pt[new_config+1] ){
		table_set();
		return;
	}
	Surface_pt[I][J]=st;
	Surface_tblX[st]=I;
	Surface_tblY[st]=J;
}
check_tree()
{
	double	new_state;
	int	i;
	
	new_state=0.0;
	for(i=64;i<128;i++) new_state+=State_tree[i];
	if( fabs(new_state-State_tree[1])/new_state > 1.0e-8 ){
		return(1);
	}

	return(0);
}
diffusion(I,J,k)
	int	I,J,k;
{
	int	new,i,j;
	
	evapolation(I,J);
	new=I_RANDOM()>>30;			/*	4 direction	*/
	I=P_X[I+x[new]+1];	/*	向前移动	*/
	J=P_Y[J+y[new]+1];
	if( k<32 ) create_A(I,J); else create_B(I,J);// 原子的类型被检查
}
choice()
{
	int	k;
	double	rnd;
	int	dmy;
123

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -