随机行走程序
INTEGER M(1:10000), NUMBER1(0:364), NUMBER2
REAL X,Y
ISEED=RTC()
DO J=1, 10000
NUMBER1=0
X=RAN(ISEED)
NUMBER1(0)=INT(365*X+1)
JJJ=1
DO I=1,364
Y=RAN(ISEED)
NUMBER2=INT(365*Y+1)
ETR=COUNT(NUMBER1.EQ.NUMBER2)
IF (ETR= =1) THEN
EXIT
ELSE
JJJ=JJJ+1
M(J)=JJJ
NUMBER1(I)=NUMBER2
END IF
END DO
END DO
DO I=1,10000
IF(M(I).LE.23) SUM=SUM+1
END DO
PRINT *,SUM/10000
END
Monte Carlo Simulation of One Dimensional Diffusion
INTEGER X,XX(1:1000,1:1000)
REAL XXM(1:1000)
! X:INSTANTANEOUS POSITION OF ATOM
! XX(J,I)：X*X ,J:第几天实验,I:第几步跳跃
! XXM(I): THE MEAN OF XX
WRITE(*,*) "实验天数JMAX ，实验次数 IMAX"
READ(*,*) JMAX,IMAX
ISEED=RTC()
DO J=1,JMAX !第几天实验
X=0 !!!
DO I=1,IMAX !第几步跳跃
RN=RAN(ISEED)
IF(RN
X=X+1
ELSE
X=X-1
END IF
XX(J,I)=X*X
END DO
END DO
OPEN(1,FILE="C:\DIF1.DAT")
DO I=1,IMAX
XXM=0.0
XXM(I)=1.0*SUM(XX(1:JMAX,I))/JMAX !!
WRITE(1,*) I, XXM(I)
END DO
! Monte Carlo Simulation of Two Dimensional Diffusion
INTEGER X,Y,XY(1:1000,1:1000)
REAL XYM(1:1000)
! X:INSTANTANEOUS POSITION OF ATOM
! XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃
! XYM(I): THE MEAN OF XY
WRITE(*,*) "实验天数JMAX ，实验次数 IMAX"
READ(*,*) JMAX,IMAX
ISEED=RTC()
DO J=1,JMAX !第几天实验
X=0 !!!
Y=0 !!!
DO I=1,IMAX !第几步跳跃
RN=RAN(ISEED)
IF(RN.LT.0.25)THEN
x=x
y=y-1
END IF
IF(RN.LT.0.5.AND.RN.GE.0.25)THEN
x=x
y=y+1
END IF
IF(RN.LT.0.75.AND.RN.GE.0.5)THEN
x=x-1
y=y
END IF
IF(RN.GE.0.75)THEN
x=x+1
y=y
END IF
XY(J,I)=X*X+Y*Y
END DO
END DO
OPEN(1,FILE="C:\DIF2.DAT")
DO I=1,IMAX
XYM=0.0
XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX !!
WRITE(1,*) I, XYM(I)
END DO
CLOSE(1)
END
! Monte Carlo Simulation of One Dimensional Diffusion
INTEGER X,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RN
REAL XYM(1:1000)
! X:INSTANTANEOUS POSITION OF ATOM
! XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃
! XYM(I): THE MEAN OF XY
WRITE(*,*) "实验天数JMAX ，实验次数 IMAX"
READ(*,*) JMAX,IMAX
XN=(/0,0,-1,1/)
YN=(/-1,1,0,0/)
ISEED=RTC()
DO J=1,JMAX !第几天实验
X=0 !!!
Y=0 !!!
DO I=1,IMAX !第几步跳跃
RN=4*RAN(ISEED)+1
X=X+YN(RN)
Y=Y+YN(RN)
XY(J,I)=X*X+Y*Y
END DO
END DO
OPEN(1,FILE="C:\DIF2.DAT")
DO I=1,IMAX
XYM=0.0
XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX !!
WRITE(1,*) I, XYM(I)
END DO
CLOSE(1)
END
做三维空间随机行走？？留作业
! Monte Carlo Simulation of One Dimensional Diffusion
INTEGER X,XY(1:1000,1:1000),y,XN(1:6),YN(1:6),ZN(1:6),RN
REAL XYM(1:1000)
! X:INSTANTANEOUS POSITION OF ATOM
! XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃
! XYM(I): THE MEAN OF XY
WRITE(*,*) "实验天数JMAX ，实验次数 IMAX"
READ(*,*) JMAX,IMAX
XN=(/0,0,-1,1,0,0/)
YN=(/-1,1,0,0,0,0/)
ZN=(/0,0,0,0,1,-1/)
ISEED=RTC()
DO J=1,JMAX !第几天实验
X=0 !!!
Y=0 !!!
Z=0
DO I=1,IMAX !第几步跳跃
RN=6*RAN(ISEED)+1
X=X+XN(RN)
Y=Y+YN(RN)
Z=Z+ZN(RN)
XY(J,I)=X*X+Y*Y+Z*Z
END DO
END DO
OPEN(1,FILE="C:\DIF2.DAT")
DO I=1,IMAX
XYM=0.0
XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX !!
WRITE(1,*) I, XYM(I)
END DO
CLOSE(1)
END
1、 随机行走程序模拟的意义
2、 演示二维随机行走一次随机抽样的结果
3、 如果随机行走是有约束的，计算方法应该如何设计
4、 讲解三维随机行走
通过该程序了解fortran 语言如何画图（通过像素画图）
USE MSFLIB
INTEGER XR ，YR ! 在的区域中画一个圆
PARAMETER XR=400,YR=400
INTEGER R, S(1:XR,1:YR)
X0=XR/2 ! 圆心位置X0，YO
Y0=YR/2
R=MIN(X0-10,Y0-10) !圆半径
S=0 !像素的初始状态（颜色）
DO I=1,XR
DO J=1,YR
IF((I-X0)**2+(J-Y0)**2
IER=SETCOLOR(S(I,J))
IER=SETPIXEL(I,J)
END DO
END DO
END
! 画一个圆（1、如何选出晶界区域；2、进一步加深对画图的理解）
USE MSFLIB
INTEGER XR ，YR ! 在的区域中画一个圆
PARAMETER XR=400,YR=400
INTEGER R, S(0:XR+1,0:YR+1), XN(1:4), YN(1:4), SNS
XN=(/0,0,-1,1/)
YN=(/-1,1,0,0/)
X0=XR/2 ! 圆心位置X0，Y0
Y0=YR/2
R=MIN(X0-10,Y0-10) !圆半径
S=0 !像素的初始状态（颜色）
DO I=1,XR
DO J=1,YR
IF((I-X0)**2+(J-Y0)**2
IER=SETCOLOR(S(I,J))
IER=SETPIXEL(I,J)
END DO
END DO
DO I=1,XR !画晶界
DO J=1,YR
NDS=0
DO K=1,4
IF(S(I,J).NE.S(I+XN(K),J+YN(K)))NDS=NDS+1
END DO
IF(NDS>0)THEN
IER=SETCOLOR(9)
ELSE
IER=SETCOLOR(8)
END IF
IER=SETPIXEL(I,J)
END DO
END DO
END
如何画有一定宽度的晶界？
！MC 模拟一个晶粒的缩小
（1、 如何定义基体和晶粒；
（2、 如何寻找边界；
（3、 如何计算能量（构造或描述问题的概率过程）
步骤：1、在基体上画一个原始晶粒，并赋予基体和原始晶粒不同的状态（取向号或体积能）
2、寻找晶界
3、MC 的一个时间步（晶粒长大过程中一种随机性）
4、计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大
① 随机选取一个初始格点；
② 若此点属于晶界，那么可以随机转变为它邻居的状态，
若不是晶界，则跳出循环，不发生转变；
③ 计算转变前后的能量变化，⊿E=⊿E v +⊿E s + ⊿E q
(自由能=体积能+表面能+能量起伏)
④ 若⊿E 小于0，则新晶相被接受，网格状态发生转变。
USE MSFLIB
PARAMETER IR=400,JR=400
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IY
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
! 定义圆心和半径
IRC=IR/2
JRC=JR/2
R=MIN(IRC,JRC)-10
! 定义基体和圆晶粒分别为状态1、状态2
IS=1
DO I=1,IR
DO J=1,JR
DISTANCE=SQRT(1.0*(I-IRC)**2+1.0*(J-JRC)**2)
IF(DISTANCE.LT.R)IS(I,J)=2
ISE=SETCOLOR(IS(I,J))
ISE=SETPIXEL(I,J)
END DO
END DO
OPEN(1,FILE="E:\LUKE.DAT")
! 寻找晶粒边界，计算能量，改变状态。
DO T=1,TMAX
DO X=1,IR
DO Y=1,JR
IX=IR*RAN(ISEED)+1
JY=JR*RAN(ISEED)+1
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1) ,IS(IX,JY+1)
,IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)
E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是圆晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
NR=8*RAN(ISEED)+1
NSTATE=ISN(NR)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTATE)
RD=RAN(ISEED)
DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25
IF(DE.LT.0.0)IS(IX,JY)=NSTATE
ISRE=SETCOLOR(IS(IX,JY))
ISRE=SETPIXEL(IX,JY)
ENDDO
ENDDO
WRITE(1,*)T,COUNT(IS.EQ.2)
ENDDO
CLOSE(1)
END
（1、 演示体积能互换的转变情况；
2、 演示IX=IR*RAN(ISEED)+1
JY=JR*RAN(ISEED)+1的替换；
3、能量关系变化；
）
作业3：基体上单晶粒形核长大
（1、 选取初始形核点
（2、 定义体积能的变化）
USE MSFLIB
PARAMETER IR=400,JR=400
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IY
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
! 定义圆心和半径
IRC=IR/2
JRC=IR/2
! 定义基体和圆晶粒分别为状态10、状态2
IS=10
IS(IRC,JRC)=2
OPEN(1,FILE="E:\LUKE.DAT")
! 寻找晶粒边界，计算能量，改变状态。
DO T=1,TMAX
DO X=1,IR
DO Y=1,JR
IX=IR*RAN(ISEED)+1
JY=JR*RAN(ISEED)+1
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)
!,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)
E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是圆晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
NR=8*RAN(ISEED)+1
NSTATE=ISN(NR)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTATE)
RD=RAN(ISEED)
DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25
IF(DE.LT.0.0)IS(IX,JY)=NSTATE
ISRE=SETCOLOR(IS(IX,JY))
ISRE=SETPIXEL(IX,JY)
ENDDO
ENDDO
WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.2))
ENDDO
CLOSE(1)
END
基体上多晶粒形核长大，
（1、 研究多晶粒长大的目的
（2、 如何定义、并区别多晶粒；
3、如何生成多晶粒形核核心；
4、为什么要有边界条件。）
步骤：1、定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等；
2、读取晶粒生长步长；
3、寻找晶界
4、计算晶界上网格点的相互作用，每个格点的相互作用
导致晶粒长大或缩小：
① 随机选取一个初始格点；
② 若此点属于晶界，那么可以随机转变为它邻居的状态，
若不是晶界，则跳出循环，不发生转变；
③ 计算转变前后的能量变化，⊿E=⊿E v +⊿E s + ⊿E q
(自由能=体积能+表面能+能量起伏)
④ 若⊿E 小于0，则新晶相被接受，网格状态发生转变。
USE MSFLIB
! 定义常数名
PARAMETER IR=400,JR=400,NMAX=100
! 定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等； INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,IY
INTEGER IGV(0:NMAX)
! 读取晶粒生长步长；
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
! 定义晶粒长大位置（100个形核点初始形核位置），并附已不同的取向号
DO I=1, NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=I
END DO
! 边界条件
IS(0,1:JMAX)=IS(IMAX,1:JMAX)
IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)
IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)
IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)
OPEN(1,FILE="E:\LUKE.DAT")
! 定义基体体积能为10, 所有晶粒体积能为1
IGV=1
IGV(0)=10
! 寻找晶粒边界。
DO T=1,TMAX
DO X=1,IR
DO Y=1,JR
IX=IR*RAN(ISEED)+1
JY=JR*RAN(ISEED)+1
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)
& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)
E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
NR=8*RAN(ISEED)+1
NSTATE=ISN(NR)
! 判断与相邻点的能量差，并决定是否改变状态
RD=RAN(ISEED)
E=COUNT(ISN.NE.NSTATE)
IG=IGV(ISN(NR))-IGV(IS(IX,JY))
DE= IG +E-E0+2.5*RD-1.25
IF(DE.LT.0.0)IS(IX,JY)=NSTATE
ISRE=SETCOLOR(MOD(IS(IX,JY),15))
ISRE=SETPIXEL(IX,JY)
ENDDO
ENDDO
WRITE(1,*)T,sqrt(1.0*COUNT(IS.NE.0))
ENDDO
CLOSE(1)
END
！iii= MOD(IS(IX,JY),15)
!If(iii==0) iii=iii+1
物理意义
1晶粒长大与时间的关系
2 形核数目与的关系
3、与EBSD 比较
4、为什么圆弧会长成近似的直线
5、柱状晶生长如何模拟呢？
Monte-Carlo 方法模拟晶粒长大过程中的新研究进展
第一个问题：采用Monte-Carlo 方法存在以下几方面的问题
（1）模拟过程只计及最近邻格点的相互作用，
难于考虑畸变场等物理场以及晶界形态对晶粒长大过程的影响，不利于在唯象层次上深入开展对晶粒
长大过程的研究。
应当指出，大部分的计算机模拟主要集中于二维系统，虽然能够在一定程度上提供晶粒长大的有关动力学和拓扑学的信息，但是，任何多晶体材料中的晶粒都是三维的，所有二维晶粒长大的模拟都是基于某种特定的假设或简化，因此模拟三维晶粒的长大更具有实际意义。同时，在目前的研究中，对晶粒生长过程的形貌演化模拟的逼真度比较好， 但在生长模型的边界处理等方面还有待完善，特别是晶粒生长与晶化温度、时间、气氛等参数密切相关，这种复杂工艺条件下的晶粒生长过程的模拟是当前该领域正待解决的难题。另外，国内外用各种方法模拟晶粒长大，主要集中在对方法本身的讨论、模拟晶粒分布函数特点及晶粒形貌等方面，而将模拟方法与实际工艺联系起来的研究者比较少。
第二个问题：算法种类：
邻居关系
1、标准算法
2、改进型算法
3、其它算法
第三个问题：模拟晶粒长大的一些应用
*柱状晶生长，树枝晶生长
1、 金属凝固过程的形核；
2、 再结晶及存在织构的再结晶（形核位置时考虑能量关系）；
3、 第二相粒子的影响（第二相粒子对晶粒长大有阻碍作用. 这种情况下模拟较复杂. 由于无法用解析解
定量描述, 用MC 方法模拟显得尤为重要；
5、异常晶粒的长大；
6、三维晶粒的长大（模拟不同的点阵结构生长过程）； 山东大学 中国有色金属学报，2008，4
7、蒙特卡洛方法不是数学的方法，而是一种多次模拟实验的方法。
元胞自动机程序（生命永不停止）
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=40000 !NMAX-随机产生的生命种子
INTEGER IS(0:1001,0:1001),IS1(0:1001,0:1001),ISN(1:8), TMAX, NUM
!IS-基体的二维数组
PRINT*,'PLEASE INPUT LOOP(TMAX)'
READ*, TMAX
ISEED=RTC()
IS=15 !"死" 的状态，基体为白色
! 赋予生命的种子，" 活" 的状态1
DO I=1, NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=1
END DO
IS1=IS
!EXECUTE THE RULE
DO T=1,TMAX
! 边界条件
IS(0,1:JMAX)=IS(IMAX,1:JMAX)
IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)
IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)
IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)
! 搜索生命存在的位置
DO IX=1,IR
DO JY=1,JR
! 判断邻居状态
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)
& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)
NUM=COUNT(ISN.EQ.1)
! 赋予生存的条件
IF((IS(IX,JY)==15.AND.NUM==3).OR.(IS(IX,JY)==1.AND
& .(NUM==3.OR.NUM== 6))) T HEN
IS1(IX,JY)=1
ELSE
IS1(IX,JY)=15
END IF
! 画图
ISRE=SETCOLOR(IS1(IX,JY))
ISRE=SETPIXEL(IX,JY)
END DO
END DO
IS=IS1
END DO
END
元胞自动机—产生和发展. 四个阶段:
1940s 诞生：Von Neumann 自我复制机.
1960-70s 起步：JH.Conway 生命游戏.
1980s 理论研究：S.Wolfram CA分类.
1980-90s 应用：HPP-FHP 格子气自动机.
C.Langton N.Packard 人工生命
元胞自动机（Cellular Automata，简称CA ，也有人称其为细胞自动机、点格自动机、分子自动机或单元自动机）是一种建立在离散的时间和空间上的动力学系统。散布在规则格网（Lattice Grid）中的每一元胞（Cell ）取有限的离散状态，遵循同样的作用规则，依据确定的局部规则作同步更新。大量元胞通过简单的相互作用而构成动态系统的演化。
与一般的动力学模型不同，元胞自动机不是由严格定义的物理方程或函数确定，而是由一系列模型构造的规则构成，凡是满足这些规则的模型都可以算做是元胞自动机模型。因此，元胞自动机是一类模型的总称，或者是一个方法框架。其特点是时间、空间、状态都是离散的，每个变量只取有限多个状态，且其状态改变的规则在时间和空间上都是局部的。 元胞自动机是一种对具有局域连通性的格点，应用局部（有时为中等范围）确定性或概率性的转换规则来描述在离散空间和时间上复杂系统演化规律的同步算法。
元胞自动机可以采取任意的维数，自动机的格点描述被认为与所选模型密切相关的基本的系统实体。独立的格点可以表示连续的体积单元，原子颗粒大原子团、汽车、人口、晶格缺陷等。
通过作用于每个格点的确定性或概率性转换规则，自动机的演化得以实施。这些规则决定了每个元胞的状态，并将其描述为其前一状态和邻近元胞的函数。在计算状态转变过程中所考虑的邻近位置的数目确定了相互作用的范围和自动机的局部演化。
元胞自动机特点：时间、空间、状态都离散的动力学模型，
利用确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化
元胞自动机模拟单晶长大
USE MSFLIB
PARAMETER IR=400,JR=400
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY !! 根据过去状态IS ，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8)
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
IRC=IR/2
JRC=JR/2
!!
! 定义基体体积能为10, 晶粒体积能为1
IS=10
IS(IRC,JRC)=1
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
! 边界条件
IS(0,1:JMAX)=IS(IMAX,1:JMAX)
IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)
IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)
IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)
!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变 DO IX=1,IR
DO JY=1,JR
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1) & ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/) E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
NR=8*RAN(ISEED)+1
NSTATE=ISN(NR)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTATE)
RD=RAN(ISEED)
IG=NSTATE-IS(IX,JY)
DE=E-E0+IG+2.5*RD-1.25
!! 用现在状态数组IS1记录边界状态的改变
IF(DE.LT.0.0)IS1(IX,JY)=NSTATE
ENDDO
END DO
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
IF(MOD(T,20).EQ.0)THEN
ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1) & ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/) ISRE=SETCOLOR(IS(IX,JY))
! 如果是边界，定义特别颜色15, 加以区分
IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0) ISRE=SETCOLOR(15) ISRE=SETPIXEL(IX,JY)
! END IF
ENDDO
ENDDO
!!
! 记录循环次数和对应的晶粒面积
WRITE(1,*)T, SQRT(1.0*COUNT(IS.EQ.1))
ENDDO
CLOSE(1)
END
元胞自动机模拟单晶长大，再附加一个规则
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=100
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY INTEGER IGV(0:NMAX)
!! 根据过去状态IS ，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8)
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
IRC=IR/2
JRC=JR/2
! 定义基体体积能为10, 晶粒体积能为1
IS=10
IS(IRC,JRC)=1
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
!!
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
!!
! 边界条件
IS(0,1:JR)=IS(IR,1:JR)
IS(IR,1:JR)=IS(iR+1,1:JR)
IS(0:IR+1,0)=IS(0:IR+1,JR)
IS(0:IR+1,JR+1)=IS(0:IR+1,1)
!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变 DO IX=1,IR
DO JY=1,JR
!!
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1) & ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/) E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
sss=0
do iii=1,8
NSTATE=ISN(iii)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTATE)
RD=RAN(ISEED)
IG=NSTATE-IS(IX,JY)
DE=E-E0+IG+2.5*RD-1.25
!! 用现在状态数组IS1记录边界状态的改变
IF(DE.LT.0.0) then
sss=sss+1
end if
end do
if(sss.gt.5) IS1(IX,JY)=11-is(ix,jy)
ENDDO
end do
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
IF(MOD(T,20).EQ.0)THEN
ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1) & ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/) ISRE=SETCOLOR(IS(IX,JY))
! 如果是边界，定义特别颜色15, 加以区分
IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0) ISRE=SETCOLOR(15) ISRE=SETPIXEL(IX,JY)
END IF
ENDDO
ENDDO
!!
! 记录循环次数和对应的晶粒面积
WRITE(1,*)T, SQRT(1.0*COUNT(IS.EQ.1))
CLOSE(1)
END
USE MSFLIB
PARAMETER IR=400,JR=400
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY !! 根据过去状态IS ，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8)
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
IRC=IR/2
JRC=JR/2
!!
! 定义基体体积能为10, 晶粒体积能为1
IS=10
IS(IRC,JRC)=1
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
! 边界条件
IS(0,1:JMAX)=IS(IMAX,1:JMAX)
IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)
IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)
IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)
!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变 DO IX=1,IR
DO JY=1,JR
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1) & ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/) num=COUNT(ISN.NE.IS(IX,JY))
! 如果不是晶粒边界，则跳出，重新循环
IF(NUM.EQ.0)CYCLE
IF((IS(IX,JY)==15.AND.NUM==3).OR.(IS(IX,JY)==1.AND
&.(NUM==2.OR.NUM== 3)).OR.NUM==8.OR.NUM==7) THEN
IS1(IX,JY)=1
ELSE
IS1(IX,JY)=15
END IF
! 随机寻找一个相邻点
! NR=8*RAN(ISEED)+1
! NSTATE=ISN(NR)
! 判断与相邻点的能量差，并决定是否改变状态
! E=COUNT(ISN.NE.NSTATE)
! RD=RAN(ISEED)
! IG=NSTATE-IS(IX,JY)
! DE=E-E0+IG+2.5*RD-1.25
!! 用现在状态数组IS1记录边界状态的改变
! IF(DE.LT.0.0)IS1(IX,JY)=NSTATE
ENDDO
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
! IF(MOD(T,20).EQ.0)THEN
ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1) & ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/) ISRE=SETCOLOR(IS(IX,JY))
! 如果是边界，定义特别颜色15, 加以区分
IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0) ISRE=SETCOLOR(15) ISRE=SETPIXEL(IX,JY)
! END IF
ENDDO
ENDDO
!!
! 记录循环次数和对应的晶粒面积
WRITE(1,*)T, SQRT(1.0*COUNT(IS.EQ.1))
ENDDO
CLOSE(1)
END
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=100
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY INTEGER IGV(0:NMAX)
!! 根据过去状态IS ，定义现在状态为IS1；为了突出边界，特别定义ISN1
INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8)
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
! 定义基体体积能为10, 晶粒体积能为1
IGV=1
IGV(0)=10
! 定义晶粒长大位置，并附能量
DO I=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IF(IS(IX0,JY0).NE.0)CYCLE
IS(IX0,JY0)=I
END DO
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
!!
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
!!
! 边界条件
IS(0,1:JMAX)=IS(IMAX,1:JMAX)
IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)
IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)
IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)
! 寻找晶粒边界，计算能量，改变状态。
!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变
DO IX=1,IR
DO JY=1,JR
!!
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)
& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)
E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
NR=8*RAN(ISEED)+1
NSTATE=ISN(NR)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTATE)
RD=RAN(ISEED)
IG=IGV(NSTATE)-IGV(IS(IX,JY))
DE=E-E0+IG+2.5*RD-1.25
!! 用现在状态数组IS1记录边界状态的改变
IF(DE.LT.0.0)IS1(IX,JY)=NSTATE
!!
!! ISRE=SETCOLOR(MOD(IS(IX,JY),15))
!! ISRE=SETPIXEL(IX,JY)
ENDDO
ENDDO
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
IF(MOD(T,20).EQ.0)THEN
ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)
& ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)
ISRE=SETCOLOR(MOD(IS(IX,JY),15))
! 如果是边界，定义特别颜色15, 加以区分
IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0) ISRE=SETCOLOR(15)
ISRE=SETPIXEL(IX,JY)
END IF
ENDDO
ENDDO
!!
! 记录循环次数和对应的晶粒面积
WRITE(1,*)T, SQRT(1.0*COUNT(IS1.gt.0)) ENDDO CLOSE(1) END
错误的多晶长大
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=100
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY
INTEGER IGV(0:NMAX)
!! 根据过去状态IS ，定义现在状态为IS1；为了突出边界，特别定义ISN1
INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8)
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
!!
! 定义基体体积能为10, 晶粒体积能为1
IGV=1
IGV(0)=10
! 定义晶粒长大位置，并附能量
DO I=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=I
END DO
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
!!
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
!!
! 边界条件
IS(0,1:JMAX)=IS(IMAX,1:JMAX)
IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)
IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)
IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)
!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变
DO IX=1,IR
DO JY=1,JR
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)
& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)
E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
sss=0
do iii=1,8
! NR=8*RAN(ISEED)+1
NSTATE=ISN(iii)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTATE)
RD=RAN(ISEED)
IG=IGV(NSTATE)-IGV(IS(IX,JY))
DE=E-E0+IG+2.5*RD-1.25
!! 用现在状态数组IS1记录边界状态的改变
! IF(DE.LT.0.0)IS1(IX,JY)=NSTATE
IF(DE.LT.0.0) then
sss=sss+1
nstate1=nstate
if(sss.gt.5) IS1(IX,JY)=nstate1
end if
end do
ENDDO
ENDDO
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
IF(MOD(T,40).EQ.0)THEN
ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)
& ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)
ISRE=SETCOLOR(MOD(IS(IX,JY),15))
! 如果是边界，定义特别颜色15, 加以区分
IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0) ISRE=SETCOLOR(15)
ISRE=SETPIXEL(IX,JY)
END IF
ENDDO
ENDDO
!!
! 记录循环次数和对应的晶粒面积
WRITE(1,*)T, COUNT(IS1.GT.0)
ENDDO
CLOSE(1)
END
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=100
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY
INTEGER IGV(0:NMAX)
!! 根据过去状态IS ，定义现在状态为IS1；为了突出边界，特别定义ISN1
INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8),ISN2(1:8)
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
!!
! 定义基体体积能为10, 晶粒体积能为1
IGV=1
IGV(0)=10
! 定义晶粒长大位置，并附能量
DO I=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=I
END DO
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
!!
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
!!
! 边界条件
IS(0,1:JMAX)=IS(IMAX,1:JMAX)
IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)
IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)
IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)
!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变
DO IX=1,IR
DO JY=1,JR
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)
& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)
E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
iiii=0
do iii=1,8
! NR=8*RAN(ISEED)+1
NSTATE=ISN(iii)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTATE)
RD=RAN(ISEED)
IG=IGV(NSTATE)-IGV(IS(IX,JY))
DE=E-E0+IG+2.5*RD-1.25
!! 用现在状态数组IS1记录边界状态的改变
! IF(DE.LT.0.0)IS1(IX,JY)=NSTATE
IF(DE.LT.0.0) then
iiii=iiii+1
isn2(iiii)=nstate
id1=iiii*ran(iseed)+1
if(iiii.gt.5) IS1(IX,JY)=isn2(id1)
else
end if
end do
ENDDO
ENDDO
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
IF(MOD(T,40).EQ.0)THEN
ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)
& ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)
ISRE=SETCOLOR(MOD(IS1(IX,JY),15))
! 如果是边界，定义特别颜色15, 加以区分
IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0) ISRE=SETCOLOR(15)
ISRE=SETPIXEL(IX,JY)
END IF
ENDDO
ENDDO
!!
! 记录循环次数和对应的晶粒面积
WRITE(1,*)T, COUNT(IS1.GT.0)
ENDDO
CLOSE(1)
END
! 相变 by Using Boundary Tracking Methord,2008,11.23
! grain growth velocity is "1.0*(IT-NT(IST))" or "a*(IT-NT(IST))**b"
! which should be not greeter than 1 cell each step.
步骤：
1、 程序的主题语句；
2、 邻居关系的不同定义方法；
3、 形核，每个时间步形核数10；
4、 获得转变点与晶粒核心的距离；
5、 比较与相场公式所得距离之差；
USE MSFLIB
PARAMETER IMAX=800,JMAX=800,NMAX=50000,NSTEP=10 ! 最多晶核数, 每步形成晶核数 INTEGER IS(0:IMAX+1,0:JMAX+1),IST(0:IMAX+1,0:JMAX+1)
INTEGER IN(1:8),JN(1:8)
INTEGER NI(1:NMAX),NJ(1:NMAX),NT(1:NMAX) !核的I,J 坐标, 形成时间
INTEGER NNS,DCN,DCNM !运行时间部，最近晶核状态, 访问CELL 与晶核的距离平方，最短距离平方 ISEED=RTC()
WRITE(*,*)"PLEASE INPUT THE TIME STEP"
READ(*,*)TMAX
IN=(/-1,-1,-1,0,0,1,1,1/)
JN=(/-1,0,1,-1,1,-1,0,1/)
OPEN(5,FILE="C:\IRX.DAT")
DO IT=1,TMAX
IS=IST
! 边界条件
IS(0,1:JMAX)=IS(IMAX,1:JMAX)
IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)
IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)
IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)
!--------------------------形核
DO I=1,NSTEP
INC=RAN(ISEED)*IMAX+1
JNC=RAN(ISEED)*JMAX+1
IF(IST(INC,JNC).NE.0)CYCLE
NN=NN+1
NI(NN)=INC
NJ(NN)=JNC
NT(NN)=IT
IST(INC,JNC)=NN !
END DO
C -------------------
DO I=1,IMAX
DO J=1,JMAX
IF(IS(I,J).GT.0)CYCLE
DCNM=2*IMAX*JMAX
NP=0
DO K=1,8
II=I+IN(K)
JJ=J+JN(K)
IF(IS(II,JJ).EQ.0)CYCLE ! 效率提高
IF(NT(IS(II,JJ)).EQ.IT)CYCLE ! 当前深刻形成的核, 当前时刻不长大
NP=NP+1 ! 只有IT 时刻以前形成的核长大
IDN=I-NI(IS(II,JJ))
JDN=J-NJ(IS(II,JJ))
IDNABS=ABS(IDN)
JDNABS=ABS(JDN)
DCN=IDN*IDN+JDN*JDN
IF(DCN.LE.DCNM)THEN
DCNM=DCN ! 访问CELL 与最近晶核的距离平方
NNS=IS(II,JJ) !
END IF
END DO
IF(NP.EQ.0)CYCLE
GR=1.0*(IT-NT(NNS)) !晶粒半径
DIS=GR-SQRT(1.0*DCNM)
IF(DIS.GE.0.0)IST(i,j)=NNS
III= MOD(IST(I,J),15)
If(III==0) III=III+1
IER=SETCOLOR(III)
IER=SETPIXEL(I,J)
END DO
END DO
! 计算相变分数，
IRX=COUNT(IST.NE.0)
WRITE(5,*)IT,1.0*IRX/(IMAX*JMAX)
END DO
CLOSE(5)
END
单晶
! 相变 by Using Boundary Tracking Methord,2008,11.13
! grain growth velocity is "1.0*(IT-NT(IST))" or "a*(IT-NT(IST))**b" ! which should be not greeter than 1 cell each step.
USE MSFLIB
PARAMETER IMAX=400,JMAX=400
INTEGER IST(0:IMAX+1,0:JMAX+1)
INTEGER DCNM ! 访问CELL 与晶核的距离平方
ISEED=RTC()
WRITE(*,*)"PLEASE INPUT THE TIME STEP"
READ(*,*)TMAX
OPEN(5,FILE="C:\IRX.DAT")
IRC=IMAX/2
JRC=JMAX/2
IST=1
IST(IRC,JRC)=2
DO IT=1,TMAX
! 边界条件
IST(0,1:JMAX)=IST(IMAX,1:JMAX)
IST(IMAX+1,1:JMAX)=IST(1,1:JMAX)
IST(0:IMAX+1,0)=IST(0:IMAX+1,JMAX)
IST(0:IMAX+1,JMAX+1)=IST(0:IMAX+1,1)
C -------------------
DO I=1,IMAX
DO J=1,JMAX
IDN=I-IRC
JDN=J-JRC
IDNABS=ABS(IDN)
JDNABS=ABS(JDN)
IF(IDNABS.GT.IMAX/2) IDN=IMAX-IDNABS !两点的最近距离 IF(JDNABS.GT.JMAX/2) JDN=JMAX-JDNABS
DCNM=IDN*IDN+JDN*JDN
GR=1.0*(IT-1) !晶粒半径
DIS=GR-SQRT(1.0*DCNM) IF(DIS.GE.0.0) IST(i,j)=2
IER=SETCOLOR(IST(I,J)) IER=SETPIXEL(I,J)
END DO
END DO
! 计算相变分数，
WRITE(5,*)IT,SQRT(1.0*COUNT(IST.EQ.2)) END DO
CLOSE(5)
END