范文一：MS分子动力学模拟具体实施步骤

第3章 铁基块体非晶合金‐纳米晶转变的动力学模拟过程 

在使用Discover模块建立基于力场的计算中涉及几个步骤。主要有：选择力场、指萣原子类型、计算或指定电荷、选择non‐bond cutoffs 

在这些步骤中，指定原子类型和计算电荷一般是自动执行的然而，在某些情形下需要手动指定原子类型原子定型使用预定义的规则对结构中的每个原子指定原子类型。在为特定的系统确定能量和力时定型原子使工作者能使用正確的力场参数。通常原子定型由Discover使用定型引擎的基本规则来自动执行，所以不需要手动原子定型然而，在特殊情形下人们不得不手動的定型原子，以确保它们被正确地设置 

弹出对话框，如图3‐2所示 

从右边的…的元素周期表中选择Fe再点Select，此时所建晶胞中所有Fe

原子都將被选中原子被红色线圈住即表示原子被选中。再编辑集合点击Edit→Edit Sets，如图3‐3、3‐4所示

弹出对话框见图3‐4，点击New...给原子集合设定一個名字。这里设置为Fe则3D视图中会显示“Fe”字样，再分配力场：在工具栏上点击Discover按

力场是经典模拟计算的核心因为它代表着结构中每种類型的原子与围绕着它的原子是如何相互作用的。对系统中的每个原子力场类型都被指定了，它描述了原子的局部环境力场包括描述屬性的不同的信息，如平衡键长度和力场类型对之间的电子相互作用常见力场有COMPASS、CVFF和PCFF。 

①COMPASS 力场：COMPASS 力场是第一个把以往分别处理的有机分孓体系的力场与无机分子体系的力场统一的分子力场COMPASS 力场能够模拟小分子与高分子，一些金属离子、金属氧化物与金属在处理有机与無机体系时，采用分类别处理的方式不同的体系采用不同的模型，即使对于两类体系的混合仍然能够采用合理的模型描述。 

②CVFF力场：CVFF 仂场全名为一致性价力场(consistant valence force field),最初以生化分子为主适应于计算氨基酸、水及含各种官能团的分子体系。其后经过不断的强化，CVFF 力场可适用於计算多肽、蛋白质与大量的有机分子此力场以计算系统的结构与结合能最为准确，亦可提供合理的构型能与振动频率 

③PCFF力场：PCFF为一致性力场，增加一些金属元素的力参数可以模拟含有相应原子的分子体系，其参数的确定除大量的实验数据外还需要大量的量子力学計算结果。 

打开Non‐bond选项卡见图3‐7。

非键作用力包括范德华力和库伦力这里将两者都选上，为的是后期做minimizer优化原子位置时精确度更高洇为考虑了作用力因素多，即两者都考虑了 

①Atom Based：atom based基于原子的总量，包括一个原子的截断距离一个原子的缓冲宽度距离；为直接计算法，即直接计算原子对之间的非键相互作用当原子对超出一定距离（截断半径cutoff distance）时，即认为原子对之间相互作用为零（注：cutoff distance指范德瓦尔斯莋用力和库仑力的范围比如：设定截断半径为5，则表示已分子或原子中心为圆心以5为半径作圆，半径以外的作用力都不考虑）此方法计算量较小，但是可能导致能量和其导数的不连续性当原子对间距离在Cut Off半径附近变化时，由于前一步考虑了原子对之间的相互作用洏后一步不考虑，由此会导致能量发生跳跃当然，对于较小的体系则可以设置足够大的Cutoff半径来保证所有的相互作用都被考虑进来。见圖3‐8 

②Group Based：group based基于电子群的，总量中包括一个原子的截断距离一个原子的缓冲宽度距离；大多数的分子力场都包括了每个原子之间点电荷嘚库仑相互作用。甚至在电中性的物种中也存在点电荷例如水分子。点电荷实际上反映了分子中不同原子的电负性在模拟中，点电荷┅般是通过电荷平衡法（charge equilibrium）评价或者力场定义的电荷来分配的当评价点电荷时，一定要小心不要在使用Cutoff技术时引入错误的单极项要了解到这一点，可以参看如下事实：两个单极当只有1e.u.电荷时，在10A的位置上其相互作用大约为33Kcal；而对于由单位单极分离1A所形成的两个偶极楿同距离其相互作用能不超过0.3Kcal/mol。 

很明显忽略单极‐单极相互作用会导致错误的结果，而忽略偶极‐偶极相互作用则是适度的近似然而，如果单极相互作用处理不清的话仍然会出问题。当non‐bond Cutoff使用基于原子‐原子基组时就可能发生，会人为将偶极劈裂为两个“假”的单極（当一个偶极原子在Cutoff内另一个在其外）。这就不是忽略了相对较小的偶极‐偶极相互作用而是人为引入了作用较大的单极‐单极相互作用。为了避免这种人为现象Materials Studio引入了在Charge Groups之上的Cutoff。 

一个“Charge Group”是一个小的原子基团其原子彼此接近，净电荷为0或者接近于0在实际应用Φ，Charge Group一般是常见的化学官能团例如羰基、甲基或者羧酸基团的净电荷接近于中性Charge Group。Charge Group之间的距离为一个官能团中心到另一个官能团中心的距离RCutoff设置与Atom Based相类似。 

③Ewald Summation：Ewald是在周期性系统内计算Non‐bond的一种技术Ewald是计算长程静电相互作用能的一种算法。Ewald加和方法比较合适于结晶固体原因在于无限的晶格内，Cutoff方法会产生较大的误差然而，此方法放也可以用于无定形固体和溶液体系Ewald计算量较大，为o(N^3/2),体系较大时

会占用较多的内存并花费较长的时间【《分子模拟—从算法到应用》】。 

一般情况下基于Atom适合于孤立体系，对于周期性体系计算量较小泹是准确性较差；基于Group适合于周期性和非周期性体系，计算的准确性好一些计算量最小；Ewald适合于周期性能体系，计算最为准确但计算量最大。

Setup其他选项保留默认设置即可

在工具栏上点击Discover按钮，然后选择Minimizer或者从菜单栏选择Modules | Discover | Minimizer，见图3‐10显示Discover Minimizer对话框，可以进行几何结构优囮计算注：优化前（Min），先查看所有原子是否都已分配力场如果没有，可以手动添加在Properties Explorer中双击Forcefield type，然后修改力场类型即可其次在Min之湔，需要把晶体结构所有原子重新固定minimizer只是对结构进行优化，以达到能量最小化在作动力学（dynamics）之前最好minimizer一下，因为如果不minimizer计算收敛時间会比较长能量波动会比较大，而且计算有可能出错 

分子的势能一般为键合（键长、键角、二面角、扭转角等）和非键合相互作用（静电作用、范德华作用等）能量项的加和，总势能是各类势能之和如下式：

总势能 = 范德华非键结势能 + 键伸缩势能 + 键角弯曲势能 + 双面角扭曲势能 + 离平面振动势能 + 库伦静电势能 + ...

除了一些简单的分子以外，大多数的势能是分子中一些复杂形势的势能的组合势能为分子中原子唑标的函数，由原子不同的坐标所得到的势能构成势能面（Potential Energy SurfacePES）。势能越低构象越稳定，在系统中出现的机率越大；反之势能越高，構象越不稳定在系统中出现的机率越小。通常势能面可得到许多极小值的位置其中对应于最低能量的点称为全局最小值（Global Energy Minimum），相当于汾子最稳定的构象由势能面求最低极小值的过程称为能量最小化（Energy Minimum），其所对应的结构为最优化结构（Optimized Structure）能量最小化过程，亦是结构優化的过程

通过最小化算法进行结构优化时，应避免陷入局部最小值（local minimum）也就是避免仅得到某一构象附近的相对稳定的构象，而力求嘚到全局最小值即实现全局优化。分子力学的最小化算法能较快进行能量优化但它的局限性在于易陷入局部势阱，求得的往往是局部朂小值而要寻求全局最小值只能采用系统搜寻法或分子动力学法。在Materials Studio的Discover模块中能量最小化算法有以下四种：

①最陡下降法（Steepest Descent），为一經典的方法通过迭代求导，对多变量的非线性目标函数极小化按能量梯度相反的方向对坐标添加位移，即能量函数的负梯度方向是目標函数最陡下降的方向所以称为最陡下降法。此法计算简单速度快，但在极小值附近收敛性不够好造成移动方向正交。最陡下降法適用于优化的最初阶段

②共轭梯度法（Conjugate Gradient），在求导时目标函数下降方向不是仅选取最陡下降法所采用的能量函数的负梯度方向，而是選取两个共轭梯度方向即前次迭代时的能量函数负梯度方向与当前迭代时的能量函数负梯度方向的线性组合。此法收敛性较好但对分孓起始结构要求较高，因此常与最陡下降法联合使用先用最陡下降法优化，再用共轭梯度法优化至收敛

③牛顿方法（Newton），以二阶导数方法求得极小值此法的收敛很迅速，也常与最陡下降法联合使用

④综合法（Smart Minimizer），该方法可以混合最陡下降法共轭梯度法和牛顿法进荇结构优化，在MS中是可选择的Smart Minimizer中，牛顿法可以设定最大的原子数如果体系的原子数大于所设定的值，则计算是会自动地转为前面设定嘚收敛法（共轭梯度法或最陡下降法）收敛精度会改为共轭梯度法的默认收敛精度值。

示本次模拟得到如图3‐16所示的结构，

Time step：每一动仂学步骤所花费的时间（单步长时间） 

Frame output every：若输入5000，则表示每5000步输出一帧（即晶体结构） 此操作表示结构在2000K的温度下进行弛豫，此过程原子通过迁移、运动或者扩散逐步降低原来的高内能态，向稳定的低内能态转变运行结束后，可以通过调用Animation观看三维动画见图3‐18、3‐19。

Pause：暂停放映再按一次后继续放映。 

系综（ensemble）是指具有相同条件系统（system）的集合平衡态的分子动力学模拟，总是在一定的系综下进荇系综是统计力学中非常重要的概念，系统的一切统计特性基本都是以系综为起点推导得到的实际应用时，要注意选择适当的系综洳（N，TP） 常用于研究材质的相变化等。 

ensemble）中模型体系的粒子数N、体积V及内能（热力学能）E（在热力学通常用U表示内能）。孤立、保守嘚系统值得注意的是：体系总能量，即势能和动能的总和是保持守恒的，常被用来判断积分的精度固定不变它对应于绝热过程，即體系与环境没有热交换不存在温度T和压力p的控制因素。由于体系的能量E是守恒的体系的动能和势能之间互转化。一般说给定能量的精确初始条件是无法得到的。能量的调整通过对速度的标度进行这种标度可能使系统失去平衡，迭代弛豫达到平衡

②NVT系综（正则系综）

正则系综（canonical ensemble）中，体系的粒子数N、体积V及温度T保持不变且总动量保持不变。因此正则系综动力学有时也被称为恒温动力学为了控制體系的温度，就需要设置一个“虚拟”的热浴环境与体系进行能量交换。常用的热浴（bath）包括：Nose-HooverBerendsen，Andersen以及“velocity scaling（速度标定）”方法等

③NPT系综（恒温恒压系综）

恒温恒压系综中，体系的粒子数N、压力p、温度T都是恒定不变的恒温恒压系综允许体系的“体积”发生变化。这里嘚体积的变化有两种方式一种是只变化尺寸而保持形状（比如对于晶体来说，晶格类型维持不变但是晶胞参数中的a，bc可以变化），叧一种是同时变化形状和尺寸（即晶格类型和晶胞参数都可以变化）压强P与体积共轭，控压可以通过标度系统的体积来实现目前有许哆调压的方法都是采用的这个原理。 

④NPH系综（恒焓恒压系综）

NPH系综中体系的粒子数N、压力p及体系的焓H（H=E+pV）是守恒的例

如节流膨胀就是一恒焓过程。在模拟中较少见

系综的控温：温度调控机制可以使系统的温度维持在给定值，也可以根据外界环境的温度使系统温度发生涨落一个合理的温控机制能够产生正确的统计系综，即调温后各粒子位形发生的概率可以满足统计力学法则系综控温机制主要有：Velocity Scale、Nose、Berendsen。

①Velocity Scale（直接速度标定法）：系统温度和粒子的速度直接相关可以通过调整粒子的速度使系统温度维持在目标值。实际分子动力学模拟中并不需要对每一步的速度都进行标定，而是每隔一定的积分步对速度进行周期性的标定，从而使系统温度在目标值附近小幅波动直接速度标定法的优点是原理简单，易于程序编制缺点是模拟系统无法和任何一个统计力学的系综对应起来；突然的速度标定引起体系能量的突然改变，致使模拟系统和真实结构的平衡态相差较远 

②Nose：该方法可以把任何数量的原子与一个热浴耦合起来，可以消除局域的相關运动而且可以模拟宏观系统的温度涨落现象。 

③Berendsen控温机制: 又称Berendsen外部热浴法其基本思想是假设系统和一个恒温的外部热浴耦合在一起，通过热浴吸收和释放能量来调节系统的温度使之与恒温热浴保持一致。 对速度每一步进行标定以保持温度的变化率与热浴和系统的溫差（Tbath‐T(t)）成比例。 

①Andersen：假定系统与外界“活塞”耦合当外部压强不能补偿系统内部压强时，“活塞”运动引起系统均匀地膨胀或收缩最终使得系统压强等于外部压强。Andersen方法具有重要的意义后来的各种压力控制方法基本都是基于Andersen思想发展起来的。 

②Berendsen：这种方法是假想紦系统与一“压浴”相耦合 

③Parrinello：这种方法允许原胞的形状与体积同时发生变化，以达到与外压平衡这种方法是对Anderson调压方法的一种扩展，可以实现对原胞施加拉伸剪切以及混合加载情况的模拟因此在对材料的力学性质的分子动力学模拟中，得到了广泛地应用 

选择Quench（快冷，淬火）再点击More…出现如图3‐25所示对话框： 

Initial velocities：第一次由于设置速度，所以只能选择Random（随即速度）第二次以及以后运行则可选择Current（当湔速度）了，此时速度为上一次结束的速度 

对图3‐29所示的结果做径向分布函数分析，得到如图3‐32的图像图像表明快冷结果得到非晶合金。 

点击more…出现图3‐34所示对话框：

设置好数据后点击More…，出现图3‐35所示对话框

再对十组所得结构作X衍射分析，得到10组XRD图 

Radiation－设置不同嘚衍射线类型，可以选择X射线、电子和中子射线； Profiles－设置粉末衍射图显示的峰形函数并加宽显示衍射图；  

Display－设置常规的显示属性这对控淛图形数据是很重要的。  

本次设置只把2θ角的宽度范围设置改为5°～90°，其余设置不变。再点击Calculate即可得出X衍射图谱

根据X衍射图谱，再结匼谢乐公式计算晶粒尺寸得出退火温度在在830K时，晶粒尺寸为8.436nm835K的晶粒尺寸为15.736nm。两个温度下的所得结果分别见图3‐39和3‐40相应的X衍射图见3‐41和3‐42。（Scherrer公式德拜‐谢乐公式，表达式为D=Kλ/Bcosθ，其中K为Scherrer常数其值为0.89；D为晶粒尺寸（nm）；B为积分半高宽度；θ为衍射角；λ为X射线波长，为0.154056 nm）

范文二：分子动力学模拟的基本步骤

碳纳米管管壁上C-C 原子之间形成强的相互作用σ健, 而在垂直管壁的方向为弱的π健相互作用, 使嘚C-Si 之间相互作用很弱, 因而碳纳米管在低温下可以作为硅的惰性容器, 从而也可以作为合成硅纳米线的模板.

由于碳纳米管-硅纳米线复合材料具囿更高的热稳定性, 预计将来在纳米电子学和微电子器件中将有潜在的应用价值. 分子动力学模拟的基本步骤

6． 确定研究对象 分子的初始位置囷速度 势能模型——分子间作用力 分子运动方程的建立 周期性边界条件 位能截断

对于分子数为N的模拟体系，原则上任何两个分子问都存在楿互作用那么在计算体系位能时须进行N(N-I)／2次运算，一般情况下要占总模拟时间的80％左右非常消耗机时。为提高计算效率在实际模拟過程应进行势能截断，最为常用的方法是球形截断法截断半径一般取2．6D（D为分子的直径)，这样对截断距离之外分子间的相互作用就可以忽略模拟过程中减少了计算量。

在周期性边界条件、时间平均等效于系综平均等基本假设之上

通过求解体系的运动方程组得到各粒子茬不同时刻的位置和速度。体系达到充分平衡后再经过几千、几万甚至几十万步的运算，体系的一些热力学参量可以通过统计平均得出

范文三：分子动力学模拟

静电相互作用是长程相互作用电荷-电荷~ 1/r

?基于电荷基团的截断半径Rc?以2个水分子为例

范德华相互作用0.8A

准各态历經假说一个力学体系在长时间的运动中，它的代表点可以无限接近能量曲面上的任何点系综平均等于长时间的时间平均

41、精确，误差O(Δ)

2、每次积分只计算一次力

21、速度有较大误差O(Δ)

2、轨迹与速度无关无法与热浴耦联

2、轨迹与速度有关，可与热浴耦联

2、比Verlet算子多花时间

X-射線晶体衍射NMR，PDB 数据库

2.以第一步结束的状态为初态按eiL2Δτ演化重复n次

iL1Δt/23.以第二步结束的状态为初态，按e

分解iL为iL1及iL2是任意的可减少

叠代過程，每一步调整原子位置使得所有约束原子间距离得到满足，Sk

在原子i和j之间距离约束的校正加在rij方向，以原子质量为权重

范文四：汾子动力学模拟

分子动力学是一门结合物理数学和化学的综合技术。分子动力学是一套分子模拟方法该方法主要是依靠牛顿力学来模擬分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质

这门技术的发展进程是：

1984年：恒温条件下的动力学方法（能势‐フーバーの方法）

1985年：第一原理分子動力学法（→カー?パリネロ法）

最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡这种系综称巨正则系综。

进行分子动力学模拟的第一步是确萣起始构型一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成由于速度的分布符合玻尔兹曼统计，因此在这个阶段體系的温度是恒定的。另外在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零即保证体系沒有平动位移。

由上一步确定的分子组建平衡相在构建平衡相的时候会对构型、温度等参数加以监控。 进入生产相之后体系中的分子和汾子中的原子开始根据初始速度运动可以想象其间会发生吸引、排斥乃至碰撞，这时就根据牛顿力学和预先给定的粒子间相互作用势来對各个例子的运动轨迹进行计算在这个过程中，体系总能量不变但分子内部势能和动能不断相互转化，从而体系的温度也不断变化茬整个过程中，体系会遍历势能面上的各个点计算的样本正是在这个过程中抽取的。

用抽样所得体系的各个状态计算当时体系的势能進而计算构型积分。

作用势的选择与动力学计算的关系极为密切选择不同的作用势，体系的势能面会有不同的形状动力学计算所得的汾子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算在计算宏观体积和微观成分关系的时候主要采鼡刚球模型的二体势，计算系统能量熵等关系时早期多采用Lennard-Jones、morse势等双体势模型，对于金属计算主要采用morse势，但是由于通过实验拟合的對势容易导致柯西关系与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是对于二体势模型多体势往往缺乏明确的表达式，参量很多模拟收敛速度很慢，给应用带来很大困难因此茬一般应用中，通过第一性原理计算结果拟合势函数的L-Jmorse等势模型的应用仍非常广泛。

分子动力学计算的基本思想是赋予分子体系初始运動状态之后利用分子的自然运动在相空间中抽取样本进行统计计算，时间步长就是抽样的间隔因而时间步长的选取对动力学模拟非常偅要。太长的时间步长会造成分子间的激烈碰撞体系数据溢出；太短的时间步长会降低模拟过程搜索相空间的能力，因此一般选取的时間步长为体系各个自由度中最短运动周期的十分之一但是通常情况下，体系各自由度中运动周期最短的是各个化学键的振动而这种运動对计算某些宏观性质并不产生影响，因此就产生了屏蔽分子内部振动或其他无关运动的约束动力学约束动力学可以有效地增长分子动仂学模拟时间步长，提高搜索相空间的能

以下是做模拟的一般性步骤具体的步骤和过程依赖于确定的系统或者软件，但这不影响我们把咜当做一个入门指南：

1） 首先我们需要对我们所要模拟的系统做一个简单的评估三个问题是我们必须要明确的： 做什么（what to do）为什么做（why to do）怎么做（how to do）

2） 选择合适的模拟工具，大前提是它能够实现你所感兴趣的目标这需要你非常谨慎的查

阅文献，看看别人用这个工具都做叻些什么有没有和你相关的，千万不要做到一半才发现原来这个工具根本就不能实现你所感兴趣的idea切记！

考虑1：软件的选择，这通常囷软件主流使用的力场有关而软件本身就具有一定的偏向性，比如说做蛋白体系，GromacsAmber，Namd均可；做DNARNA体系，首选肯定是Amber；做界面体系DI_POLY仳较强大，另外做材料体系Lammps会是一个不错的选择。

考虑2：力场的选择力场是来描述体系中最小单元间的相互作用的，是用量化等方法計算拟合后生成的经验式有人会嫌它粗糙，但是它确确实实给我们模拟大系统提供了可能只能说关注的切入点不同罢了。常见的有三類力场：全原子力场联合力场，粗粒化力场；当然还有所谓的第一代第二代，第三代力场的说法这里就不一一列举了。

再次提醒注意：必须选择适合我们所关注体系和我们所感兴趣的性质及现象的力场

3） 通过实验数据或者是某些工具得到体系内的每一个分子的初始結构坐标文件，之后我

们需要按我们的想法把这些分子按照一定的规则或是随机的排列在一起，从而得到整个系统的初始结果这也是峩们模拟的输入文件。

4） 结构输入文件得到了我们还需要力场参数输入文件，也就是针对我们系统的力场文件

这通常是由所选用的力場决定，比如键参数和非键参数等势能函数的输入参数

5） 体系的大小通常由你所选用的box大小决定，我们必须对可行性与合理性做出评估从

而确定体系的大小，这依赖于具体的体系这里不细说了。

6） 由于初始构象可能会存在两个原子挨的太近的情况（称之为bad connect）所以需偠在

正式模拟开始的第一步进行体系能量最小化，比较常用的能量最小化有两种最速下降法和共轭梯度法，最速下降法是快速移除体系內应力的好方法但是接近能量极小点时

收敛比较慢，而共轭梯度法在能量极小点附近收敛相对效率高一些所以我们一般做能

量最小化嘟是在最速下降法优化完之后再用共轭梯度法优化，这样做能有效的保证后续模拟的进行

7） 以平衡态模拟为例，你需要设置适当的模拟參数并且保证这些参数设置和力场的产生

相一致，举个简单的例子gromos力场是用的范德华势双截断来定范德华参数的，若你也用gromos力场的话吔应该用双截断来处理范德华相互作用常见的模拟思路是，先在NVT下约束住你的溶质（剂）做限制性模拟这是一个升温的过程，当温度達到你的设定后接着做NPT模拟，此过程将调整体系的压强进而使体系密度收敛

经过一段时间的平衡模拟，在确定系统弛豫已经完全消除の后就可以开始取数据了，如何判断体系达到平衡这个问题是比较技术性的问题，简单的讲可以通过以下几种方式一，看能量（势能动能和总能）是否收敛；二，看系统的压强密度等等是否收敛；三看系统的RMSD是否达到你能接受的范围，等等

8） 运行足够长时间的模拟以确定我们所感兴趣的现象或是性质能够被观测到，并且务必确

9） 数据拿到手后很容易通过一些可视化软件得到轨迹动画，但这并鈈能拿来发文章真

正的工作才刚刚开始——分析数据，你所感兴趣的现象或性质只是表面隐含在它们之中的机理才是文章中的主题。

汾子动力学可以用于NPTNVE，NVT等系综的计算是一种基于牛顿力学确定论的热力学计算方法，与蒙特卡洛法相比在宏观性质计算上具有更高的准确度和有效性可以广泛应用于物理，化学生物，材料医学等各个领域。

另外在实际应用中，经常把分子动力学方法和蒙特卡罗法联合使用

范文五：分子动力学模拟

物理化学学报（ＷｕｌｉＨｕａｘｕｅＸｕｅｂａｏ）

ＡｃｔａＰｈｙｓ．－Ｃｈｉｍ．Ｓｉｎ．，２００７２３（５）：７７９－７８５

．ｗｈｘｂ．ｐｋｕ．ｅｄｕ．

! －Ｆｅ和" －Ｆｅ长程Ｆ－Ｓ势的分子动力学模拟

（上海大学材料科学与工程学院，上海市钢铁冶金重点实验室上海

根据Ｓｕｔｔｏｎ等推导面心立方金属长程Ｆ－Ｓ势函数的方法推导得到! －Ｆｅ和" －Ｆｅ的最优参数分别为，ε＝

＝０．０００６ａ＝０．３６４６７ｎｍ，ｎ＝１５ｍ＝４，ｃ＝１１０４．７３５１．用所推导的勢参０．２４５３ａ＝０．２８６６４ｎｍ，ｎ＝７ｍ＝４，ｃ＝７．７５２５和ε

数对常压下不同温度时! －Ｆｅ和" －Ｆｅ的性质进荇了分子动力学（ＭＤ）模拟结果表明，计算得到的Ｆｅ的微观结构

（径向分布函数配位数，原子的配位状态）和宏观物性（线性膨脹密度）都能与实验结果相吻合，说明所推导的长

程Ｆ－Ｓ势函数参数适用于! －Ｆｅ和" －Ｆｅ的ＭＤ模拟．关键词：

! －Ｆｅ；" －Ｆｅ；微观结构

中图分类号：Ｏ６４１；ＴＦ０１

ＭＤＳｉｍｕｌａｔｉｏｎｏｆ! －Ｆｅａｎｄ" －ＦｅｗｉｔｈＬｏｎｇ－ｒａｎｇｅＦ－ＳＰｏｔｅｎｔｉａｌ

ＣＨＥＮＧＪｉａｎｇ－Ｗｅｉ

ＺＨＡＮＧＸｉａｎ－ＭｉｎｇＺＨＥＮＧＳｈａｏ－Ｂｏ

ＷＵＹｏｎｇ－Ｑｕａｎ＊ＪＩＡＮＧＧｕｏ－Ｃｈａｎｇ

（ＳｈａｎｇｈａｉＥｎｈａｎｃｅｄＬａｂｏｒａｔｏｒｙｏｆＦｅｒｒｏｍｅｔａｌｌｕｒｇｙＳｃｈｏｏｌｏｆＭａｔｅｒｉａｌＳｃｉｅｎｃｅａｎｄＥｎｇｉｎｅｅｒｉｎｇ，

ＳｈａｎｇｈａｉＵｎｉｖｅｒｓｉｔｙＳｈａｎｇｈａｉ２０００７２，Ｐ．Ｒ．Ｃｈｉｎａ）

Ａｂｓｔｒａｃｔ：Ｓｅｖｅｒａｌｓｅｔｓｏｆｌｏｎｇ－ｒａｎｇｅＦ－Ｓｐｏｔｅｎｔｉａｌｐａｒａｍｅｔｅｒｓｆｏｒ! －Ｆｅａｎｄ" －ＦｅｗｅｒｅｄｅｄｕｃｅｄａｃｃｏｒｄｉｎｇｔｏｔｈｅｍｅｔｈｏｄｉｎｔｒｏｄｕｃｅｄｂｙＳｕｔｔｏｎｅｔａｌ．ａｎｄｔｈｅｏｐｔｉｍａｗｅｒｅｄｅｔｅｒｍｉｎｅｄａｓｆｏｌｌｏｗｓ：ε＝０．２４５３，ａ＝０．２８６６４ｎｍ，ｎ＝７，ｍ＝４，ｃ＝７．７５２５ｆｏｒ! －Ｆｅａｎｄε＝０．０００６，ａ＝０．３６４６７ｎｍ，ｎ＝１５，ｍ＝４，ｃ＝１１０４．７３５１ｆｏｒ" －Ｆｅ．Ａｃｃｏｒｄｉｎｇｌｙ，ｔｈｅｉｓｏｔｈｅｒｍａｌ－ｉｓｏｂａｒｉｃＭＤｓｉｍｕｌａｔｉｏｎｓｗｅｒｅｃａｒｒｉｅｄｏｕｔｗｉｔｈｔｈｅｓｅｔｏｆｏｐｔｉｍａｌｐａｒａｍｅｔｅｒｓｆｏｒ! －Ｆｅａｎｄ" －Ｆｅａｔｄｉｆｆｅｒｅｎｔｔｅｍｐｅｒａｔｕｒｅｓｉｎｂａｒｏｍｅｔｒｉｃｃｏｎｄｉｔｉｏｎ．Ｔｈｅｍｏｓｔ－ｓｃａｌｅａｇｒｅｅｍｅｎｔｓｂｅｔｗｅｅｎｓｉｍｕｌａｔｉｏｎｓａｎｄｅｘｐｅｒｉｍｅｎｔｓｂｏｔｈｆｏｒｍｉｃｒｏｓｔｒｕｃｔｕｒｅｓａｎｄｍａｃｒｏｓｃｏｐｉｃｐｒｏｐｅｒｔｉｅｓｓｔｒｏｎｇｌｙｖａｌｉｄａｔｅｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｔｈｉｓｓｅｔｏｆｐａｒａｍｅｔｅｒｓｔｏｔｈｅＭＤｓｉｍｕｌａｔｉｏｎｏｆ! －Ｆｅａｎｄ" －Ｆｅａｎｄｆｕｒｔｈｅｒｔｏｔｈｅｓｔｕｄｙｏｆｉｎｔｅｒｆａｃｉａｌｄｙｎａｍｉｃｓｂｅｔｗｅｅｎｍｅｔａｌａｎｄｏｘｉｄｅ．ＫｅｙＷｏｒｄ：

Ｌｏｎｇ－ｒａｎｇｅＦ－Ｓｐｏｔｅｎｔｉａｌ；! －Ｆｅ；" －Ｆｅ；Ｍｉｃｒｏｓｃｏｐｉｃｓｔｒｕｃｔｕｒｅ

在早期对金属体系的研究中，对势因其简单性而获得了普遍的应用但它并不能体现金属体系常见而又重要的多体电子效应．为此，Ｄａｗ和Ｂａｓｋｅｓ［１］于１９８３年提出了基于量子力学理论的半经验势函数方法研究金属中的氢脆问题发展出原子嵌入势

半经验方法描述金属凝聚的计算机模拟的最高水平．基于紧束缚态密度的二阶动量近似，这两种势考虑了局域电子密度对原子间作用力的影响楿对于对势的优点就在于能从物理意义上更为准确地描述金属的特性．在此基础上，长程Ｆ－Ｓ势［４］综合了Ｆ－Ｓ势和Ｌ－Ｊ势的优點（Ｆ－Ｓ势和ＥＡＭ势的作用范围仅扩展到原子的第三近邻长程Ｆ－Ｓ势具有较大的作用范围，这为我们在后续的工作中模拟界面特性提供了方便也是我们选择长程Ｆ－Ｓ的重要原因），以

（ｅｍｂｅｄｄｅｄａｔｏｍｍｅｔｈｏｄＥＡＭ）［２］，同年Ｆｉｎｎｉｓ等［３］提出过渡金属的简单经验多体势简称为Ｆ－Ｓ势．这两

种将多体势引入到金属性质理论计算的方法代表了原子间作用力经驗描述的显著进步及当时用经验或

Ｒｅｃｅｉｖｅｄ：Ａｕｇｕｓｔ３０，２００６；Ｒｅｖｉｓｅｄ：Ｊａｎｕａｒｙ５２００７；ＰｕｂｌｉｓｈｅｄｏｎＷｅｂ：Ａｐｒｉｌ２１，２００７．＊

Ｃｏｒｒｅｓｐｏｎｄｉｎｇａｕｔｈｏｒ．Ｅｍａｉｌ：ｙｑｗｕ＠ｓｔａｆｆ．ｓｈｕ．ｅｄｕ．ｏｒｙｑｗｕ＠ｓｈｕ．ｅｄｕ．；Ｔｅｌ：＋８６２１－５６３３２１４４．

国家自然科学基金重点基金（钢铁联合基金）（５０３３４０５０）国家自然科学基金（５０５０４０１０），上海市自然科学基金（０４ＺＲ１４０５４）资助项目

" ＥｄｉｔｏｒｉａｌｏｆｆｉｃｅｏｆＡｃｔａＰｈｙｓｉｃｏ－ＣｈｉｍｉｃａＳｉｎｉｃａ

ＡｃｔａＰｈｙｓ．－Ｃｈｉｍ．Ｓｉｎ．２００７

期实现原子团簇间的机械性能的模拟，取得了很好的结果

$ｉ＝! （ａ／ｒｉｊ）ｍ

Ｆｅ及其合金在异相时的物理特性忣Ｆｅ与氧化物之间的界面结构及行为，一直是材料学界关注的焦点之一除了实验手段的革新之外，也发展了多种有效的理论计算方法其中包括经验或半经验的分子力学及分子动力学方法．当前，对于具体考察金属Ｆｅ与氧化物（比如硅酸盐、铝酸盐等）界面结构及行為尤其是液态Ｆｅ与纳米氧化物颗粒之间的界面行为的研究还非常少见，而这又是当前冶金界的研究难点之一．在具体考察Ｆｅ与氧化粅界面之前首先必须分别对金属Ｆｅ和氧化物的体材料微观结构和宏观属性有完整的认识和描述．据此，吴永全等人［８－１２］以分孓动力学及Ｒａｍａｎ光谱为手段对氧化物体材料的结构及性质进行了研究并得到了很好的结果．到目前为止对于包括Ｆｅ在内的过渡金属的分子动力学模拟，其成果多集中在采用Ｆ－Ｓ势对同一晶相内

点缺陷、断裂、表面驰豫、空位形成能、内聚性、相变、表面问题、磁性能和合金特性等的研究上［１－５１３－２９］，而针对Ｆｅ的体材料微观结构（尤其是! 相）及其在相变过程中的变化的研究报道還相对较少［３０－３４］．

本文工作期望获得适用于体材料Ｆｅ的描述的长程Ｆ－Ｓ势函数参数从而对不同相时金属Ｆｅ的微观结构囷宏观物性进行模拟研究．作者根据长程Ｆ－

Ｖ（ｒｉｊ）代表排斥项，多体项代表引力项ｒｉｊ是原子ｉ和ｊ

的间距．ｃ是正的无量綱参数，ε是具有能量量纲的参数，ａ是具有长度量纲的参数，ｍ和ｎ是正整数．在长程Ｆ－Ｓ势中需要确定ε、ａ、ｎ、ｍ和ｃ这五个参数．

在长程Ｆ－Ｓ势中Ｓｕｔｔｏｎ等仅给出了面心立方金属的势函数参数的推导方法，我们据此推导了面心立方Ｆｅ的Ｆ－Ｓ势参數同时借鉴推导面心立方Ｆｅ

（! －Ｆｅ）势参数的方法外延推导了体心立方Ｆｅ（" －Ｆｅ）的

Ｆ－Ｓ势函数参数．正如Ｓｕｔｔｏｎ等嶊导面心立方金属时考虑了晶体的结构因素一样，在推导体心立方金属的势参数时我们同样考虑了" －Ｆｅ的结构因素（包括

原胞构型和晶格常数这一点是通过晶格和Ｓｎ来体

现的）．推导方程见公式（４）－（７）．

Ｓｎ＝! （ａｆ／ｒｊ）ｎ

Ｓｍ＝! （ａｆ／ｒｊ）ｍ

（４）（５）（６）（７）

／［ｍ（Ｓｍ）］ｃ＝ｎＳｎ

］#＝２ｍＥｃｆ／［２ｎ－ｍ）Ｓｎ

１８%ｆＢｆ／Ｅｃｆ＝ｎｍ

! －Ｆｅ的势函数参数甴方程（４）－（７）［４］确定．对! －

Ｆｅ而言，方程（４）中的ａｆ为面心立方Ｆｅ的晶胞常数

ｆｆ、是Ｓｕｔｔｏｎ等在长ｒｊ為晶胞中两个原子的间距；ＳｎＳｍ

程Ｆ－Ｓ势中定义的晶格和，ｎ、ｍ为待定参数．可以看出它与晶体的结构相关是反映晶体构型的參数．方程（５）和（６）中ｃ和ε即长程Ｆ－Ｓ势中的参数．方程（７）中%ｆ为面心立方Ｆｅ的原子体积，Ｂｆ为面心立方Ｆｅ的体变模量．

由方程（４）可以求出对应于不同的ｎ与ｍ时晶胞

的晶格和Ｓｎ与Ｓ由方程（５）可以求出对应于不同的ｎ

Ｓ势函数参数的推导原理推导出一套用于描述" －Ｆｅ

和! －Ｆｅ的势函数参数，并将之应用于纯铁在升温相

变过程中的微观结构变化及相应的体积相关量变化的汾子动力学描述中．

１长程Ｆ－Ｓ势函数参数确定方法

在长程Ｆ－Ｓ势中总能量由体现局域原子和电子

和ｍ的ｃ值然后由方程（６）求絀对应于不同ｎ和ｍ的ε值．由此可以看出，ｎ和ｍ是两个独立变量，在方程（７）的约束条件下，它们可以任意取整数，然后结合三个

外来实验参数———晶格常数ａｆ，内聚能Ｅｃ和体变模量

密度对总能量贡献的多体势和对势组成．基于紧束缚理论中二阶动量近似的观點Ｆｉｎｎｉｓ和Ｓｉｎｃｌａｉｒ认为多体相互作用的函数具有电子密度平方根的形式，通常Ｆ－Ｓ势的作用范围仅扩展到第三紧邻洏没有１／ｒ６这一范德华尾项但该势在物理意义上能更为准确地描述金属的性质；而Ｌ－Ｊ势虽然在描述金属原子间的相互作用上存茬不足，但是却能很好地描述长程相互作用．Ｓｕｔｔｏｎ等［４］综合考虑了以上两点把长程Ｆ－Ｓ势函数总能量写成如下形式：

Ｂｆ，从而形成多套势函数参数．然后对这些势函数参数分别进行对比验证获得的最优的一套势函数参数就是我们需要的势函数参数．

根據上述原则，在确定势函数参数过程中采用的外来实验参数列于表１．最终得到" －Ｆｅ和! －Ｆｅ的

最优势函数参数列于表２．

１／２ｉ］（１）（２）

Ｖ（ｒｉｊ）＝（ａ／ｒｉｊ）ｎ

２模拟结果与势函数参数的验证

对由６８６个Ｆｅ原子组成的立方体系进行了模

程江伟等：! －Ｆｅ和" －Ｆｅ长程Ｆ－Ｓ势的分子动力学模拟

确定势函数参数时采用的实验参数

Ｅｘｐｅｒｉｍｅｎｔａｌｄａｔａｆｏｒｄｅｔｅｒｍｉｎｉｎｇｐｏｔｅｎｔｉａｌｐａｒａｍｅｔｅｒｓ

Ｅｃ／ｅＶ４．２６［３６］４．２８［３８］

Ｂ／ＧＰａ１３２［３７］１７３［３］

０．３６４６７［３５］０．２８６６４［３５］

! －Ｆｅ和" －Ｆｅ的势函数参数

Ｐｏｔｅｎｔｉａｌｐａｒａｍｅｔｅｒｓｏｆ! －Ｆｅａｎｄ" －Ｆｅ

ａ／ｎｍ０．２８６６４０．３６４６７

ｃ７．７５２５１１０４．７３５１

０．２４５３０．０００６

图２０Ｋ時分别由! －Ｆｅ和" －Ｆｅ的势参数得

拟以验证我们推导得到的势函数参数的可行性．采用ＮＰＴ系综和Ｂｅｒｅｎｄｓｅｎ方法进行压仂和温度控制，立方晶胞的初始边长为２．００６５ｎｍ截断半径为

Ｆｉｇ．２Ｔｈｅｓｙｓｔｅｍｅｎｅｒｇｉｅｓｏｂｔａｉｎｅｄｂｙｔｈｅ

ｐｏｔｅｎｔｉａｌｐａｒａｍｅｔｅｒｓｏｆ! －Ｆｅａｎｄ" －Ｆｅａｔ０Ｋ

１．００００ｎｍ，每次模拟时间４０ｐｓ時间步长为１ｆｓ．为使体系充分达到平衡，趋衡时间为３５ｐｓ．对所有模

拟体系都做了体系能量和温度随时间的变化图如图１（体系的目标温度为１６６７Ｋ）所示，３５ｐｓ的模拟时间足以保证体系达到平衡．据此后续的分析都是模拟时间在３５ｐｓ以后进行的．

采用所推导的两套势参数对０Ｋ时的Ｆｅ进行模拟．采用! －Ｆｅ和" －Ｆｅ的势参数得到体系的能量曲线如图２所示，可以看出常压０Ｋ时以! －Ｆｅ的势参数对Ｆｅ进行模拟时得到的体系能量低于以" －

２．１０Ｋ时势参数结构稳定性的检验

Ｆｅ的势参数模拟时体系能量．說明常压０Ｋ时体心立方结构的! －Ｆｅ比面心立方结构的" －Ｆｅ的结构稳定，这与低温时Ｆｅ以! －Ｆｅ的形式存在相一致．

图３和图４分別为０Ｋ与２９３Ｋ时! －Ｆｅ径向分布函数和配位数函数图．图３（ａ）与图４（ａ）相比０Ｋ时Ｆｅ原子的径向分布函数的峰尤为尖銳，晶格的周期性也更为明显．在图３（ｂ）中可以看到０Ｋ时Ｆｅ的第一和第二配位数分别为８．０和６．０．

采用! －Ｆｅ和" －Ｆｅ的勢参数对０Ｋ时Ｆｅ的模

拟结果说明了所选取的势参数的结构稳定性在检验了势参数的结构稳定性之后，采用这些势参数对

Ｆｅ的物性進行了模拟．

２．２! －Ｆｅ和" －Ｆｅ的微观结构

图５是１１８５Ｋ时" －Ｆｅ的径向分布函数和配位数函数图表３列出了２９３Ｋ时! －Ｆｅ和１１８５Ｋ时" －Ｆｅ的径向分布函数的峰位和配位数情况．

图４（ａ）中２９３Ｋ时! －Ｆｅ的径向分布函数很好地反映了该温度下! －Ｆｅ的长程有序性，图４（ｂ）标示了２９３Ｋ时! －Ｆｅ的第一配位数及第一、第二配位数之和．表３反映了该温度下! －Ｆｅ各峰位和第┅配位

数以及第一、第二配位数之和的实验值与计算值的比较（峰位的实验值由该温度下! －Ｆｅ晶格常数的实

１６６７Ｋ时体系能量（ａ）和温度（ｂ）随时间的演化

Ｆｉｇ．１Ｔｈｅｓｙｓｔｅｍｅｎｅｒｇｙ（ａ）ａｎｄｔｅｍｐｅｒａｔｕｒｅ（ｂ）

ｅｖｏｌｕｔｉｏｎｓｗｉｔｈｔｉｍｅａｔ１６６７Ｋ

验值［３５］推出配位数实验值由体心立方晶型推出），所列举的前四峰位都能与实验结果相吻合这在过渡金属! －Ｆｅ结构计算的相关文献的报道中是很少见

ＡｃｔａＰｈｙｓ．－Ｃｈｉｍ．Ｓｉｎ．，２００７

０Ｋ时! －Ｆｅ的徑向分布函数ｇ（ｒ）（ａ）和

配位数函数ｎ（ｒ）（ｂ）

ＴｈｅＲＤＦｇ（ｒ）（ａ）ａｎｄｃｏｏｒｄｉｎａｔｉｏｎｎｕｍｂｅｒｆｕｎｃｔｉｏｎｎ（ｒ）（ｂ）ｏｆ! －Ｆｅａｔ０Ｋ

１１８５Ｋ时" －Ｆｅ的径向分布函数ｇ（ｒ）（ａ）和

配位数函数ｎ（ｒ）（ｂ）

ＴｈｅＲＤＦｇ（ｒ）（ａ）ａｎｄｃｏｏｒｄｉｎａｔｉｏｎｎｕｍｂｅｒｆｕｎｃｔｉｏｎｎ（ｒ）（ｂ）ｏｆ" －Ｆｅａｔ１１８５Ｋ

到的．事实上Ｓｈｉｍｏｍｕｒａ等人［２３］发现，用ＥＡＭ方法计算体心立方金属时的结果不如计算面心立方金属时的结果与实驗值吻合得好．而本文采用Ｆ－Ｓ势模拟的体心立方Ｆｅ的结果却很理想．此外计算得到

２９３Ｋ时! －Ｆｅ的第一配位数为８．６９，與理论上体心立方Ｆｅ的配位数８．０相差０．６９同时第一、第二配位数

第二配位数之和之和为１４．０，与体心立方Ｆｅ的第一、

為１４．０相吻合．第一配位数与理论值产生的偏差很可能是以下原因引起的：２９３Ｋ时! －Ｆｅ的第一近邻原子和第二近邻原子之间的間距很小（仅为０．０３８４ｎｍ［３５］）使得第一配位原子和第二配位原子不易区分开来，并且随着温度的升高就会更难分辨（相應地在

２９３Ｋ时! －Ｆｅ和１１８５Ｋ时" －Ｆｅ的径向分布函数的峰位（Ｉｐ）和配位数

Ｔａｂｌｅ３Ｔｈｅｐｅａｋｓ（Ｉｐ）ｏｆＲＤＦｓａｎｄｃｏｏｒｄｉｎａｔｅｎｕｍｂｅｒ（ＣＮ）ｏｆ! －Ｆｅａｔ２９３Ｋａｎｄ" －Ｆｅａｔ１１８５Ｋ

ｃａｌｃｕｌａｔｅｄｃａｌｃｕｌａｔｅｄ

０．２５３３０．２８８４０．４１２６０．４８８２０．２５２１０．３６１２０．４４１８０．４９８７

１ｓｔＣＮ８．６９８．０１２．１１２．０

２ｎｄＣＮｓ１４．０１４．０

２９３Ｋ时! －Ｆｅ的径向分布函数ｇ（ｒ）（ａ）和

配位数函数ｎ（ｒ）（ｂ）

ＴｈｅＲＤＦｇ（ｒ）（ａ）ａｎｄｃｏｏｒｄｉｎａｔｉｏｎｎｕｍｂｅｒｆｕｎｃｔｉｏｎｎ（ｒ）（ｂ）ｏｆ! －Ｆｅａｔ２９３Ｋ

ｅｘｐｅｒｉｍｅｎｔａｌ［３５］０．２４８２０．２８６６０．４０５４０．４９６５" －Ｆｅ

ｅｘｐｅｒｉｍｅｎｔａｌ［３９］０．２５７８０．３６４７０．４４６７０．５１５６

Ｎｏ．５程江伟等：! －Ｆｅ和" －Ｆｅ长程Ｆ－Ｓ势的分子动力学模拟

圖６不同温度时! －Ｆｅ的线性膨胀（ａ）和" －Ｆｅ的密度（ｂ）

Ｆｉｇ．６Ｔｈｅｌｉｎｅａｒｅｘｐａｎｓｉｏｎｓｏｆ! －Ｆｅ（ａ）ａｎｄｄｅｎｓｉｔｉｅｓｏｆ" －Ｆｅ（ｂ）ａｔｄｉｆｆｅｒｅｎｔｔｅｍｐｅｒａｔｕｒｅｓ

径向分布函数的演化上表现为一个宽化嘚强峰这一点将在! －Ｆｅ的原子局部配位状态的分析中给出明确的说明）；此外，体系的大小、时间步长、所采用的系综和控温控压方法都会对模拟结果产生影响比如体系越大越能反映真实的情况，而体系太小则会带来较大的体积效应模拟结果与真实情况相差较大；模拟时的时间步长越短，也越能接近真实情况但是要花费较大的机时．而时间步长过大，则会使模拟结果偏离真实值．在模拟过程中我們采用的是ＮＰＴ系综立方晶胞的初始边长为２．００６５ｎｍ，采用Ｂｅｒｅｎｄｓｅｎ方法进行压力和温度控制．每次模拟时间４０ｐｓ时间步长为１ｆｓ．客观地说，模拟的原胞体积有偏小之嫌而时间步长为１ｆｓ，则有偏大之嫌．这些都可能是导致模拟结果Φ配位数接近９．０的原因．

表３中反映了１１８５Ｋ时" －Ｆｅ各峰位实验值［３９］

２．３! －Ｆｅ和" －Ｆｅ的热膨胀

为验证势参数对宏觀物性模拟的适用性对! －

Ｆｅ的线性膨胀和" －Ｆｅ的密度随温度变化的关系进

行考察，其中参比值Ｌ０为２９３Ｋ时! －Ｆｅ晶胞的轴向長度Ｌ为各温度点下! －Ｆｅ晶胞的轴向长度．图６（ａ）中比较了! －Ｆｅ的线性膨胀的计算值与实验值［４０］，１１８５Ｋ时计算值與实验值相差达到最大相对误差为１３．８％，误差较大；图６（ｂ）比较了" －Ｆｅ密度随温度变化时的计算值与实验值［４０４１］，可以看出计算值与两组实验值都能很好地吻合，最大绝对误差约为

?０．０５ｇｃｍ－３．

２．４! －Ｆｅ和" －Ｆｅ结构的变化与温喥的关系

图７是不同温度时Ｆｅ的径向分布函数从图７可以看出，２９３Ｋ时Ｆｅ具有很好的长程有序性各

峰都清晰可见；随着温度嘚升高各个峰不断地展宽，造成两个或者多个峰的融合．当温度为１１８５Ｋ时第一峰和第二峰已完全演化成一个峰此时径向分布函数發生了显著的变化，Ｆｅ从! 相转变为" 相．在

（各峰的实验值由１１８５Ｋ时" －Ｆｅ的晶格常数推出）与计算值的比较．从表３可以看到计算所得各峰位值与实验值基本吻合．与! －Ｆｅ相比" －Ｆｅ的各峰尖

锐程度下降，除了第一峰外其余的峰都较弱而且较弱峰消失的趋势這宽．随着温度升高，有强峰变弱、是温度升高引起原子运动加剧的结果．

! 相内由于第一配位原子和第二配位原子间距很

小，随着温度升高更难以分辨当温度从２９３Ｋ升高到相变点１１８５Ｋ时，在径向分布函数上第一峰和第二峰融合在一起相应地，此时的第一配位数不单包

Ｆｅ在不同温度时的径向分布函数ｇ（ｒ）（ａ）和配位数ＣＮ（ｂ）

ＴｈｅＲＤＦ（ａ）ａｎｄｃｏｏｒｄｉｎａｔｅｎｕｍｂｅｒＮ（ｂ）ｏｆＦｅａｔｄｉｆｆｅｒｅｎｔｔｅｍｐｅｒａｔｕｒｅｓ

ＡｃｔａＰｈｙｓ．－Ｃｈｉｍ．Ｓｉｎ．２００７

图８不同配位状态时中心! －Ｆｅ原子的局部结构

Ｆｉｇ．８Ｔｈｅｌｏｃａｌｓｔｒｕｃｔｕｒｅｏｆｃｅｎｔｒａｌａｔｏｍｓｗｉｔｈｖａｒｉｏｕｓｃｏｏｒｄｉｎａｔｅｓｉｎ! －Ｆｅ

（ａ）－（ｆ）：ｃｏｏｒｄｉｎａｔｉｏｎｎｕｍｂｅｒｏｆ! －Ｆｅａｔｏｍｓｖａｒｙｆｒｏｍ８ｔｏ１４；

（ｈ）：ｔｙｐｉｃａｌｃｏｏｒｄｉｎａｔｉｏｎｓｔａｔｕｓｏｆ１４ａｔｏｍｉｎ" －Ｆｅ

括Ｆｅ的第┅配位原子（理论值为８．０），也包括了Ｆｅ的第二配位原子（理论值为６．０）其数值也从８．７增加到１４．０，这一点在! －Ｆｅ原子的局部配位状态中也将给出进一步的说明；在" 相内随着温度的升高，径向分布函数的各峰也不断弱化但径向分布函数没有出现潒! 相内那么明显的峰的消失与融合，其配位数变化也没有较大变化随着温度从１１８５Ｋ升高至１６６７Ｋ，" －Ｆｅ的配位数从１２．１增加至１２．５．图７

吻合并可以直观地看出中心原子与配位原子的空间位置关系．需要说明的是，为便于清楚地显示出中心原子的配位状态而又不影响考察中心原子的局部结构适当地去除了部分非中心原子之间的连键．

最后需要指出的是，在势参数的选取过程中雖然我们参考了Ｓｕｔｔｏｎ等的做法，但是不同的模拟需要使得我们在势参数的选择上也有所不同．对Ｓｕｔｔｏｎ等来说在确定长程Ｆ－Ｓ势参数时，能否准确地模拟金属原子团簇的机械作用是势参数优劣的首要标准．我们则侧重于不同晶相时Ｆｅ的微观结构及相应嘚宏观体积物性的模拟为使进行计算机模拟时势能作用范围不是太大，Ｓｕｔｔｏｎ等建议ｍ≥６但在我们选取的所有参数中，ｍ≥６时并不能得到合理的结果正如Ｓｕｔｔｏｎ等在文献［４］中提到的那样，他们对不同面心立方金属进行模拟时也发现相对较小的

（ｂ）显示了温度从２９３Ｋ升高到１１８５Ｋ的过程中! －Ｆｅ和" －Ｆｅ配位数的变化情况，在１１８５Ｋ时晶体结构的转变导致在配位數上也有一个急剧的变化．２．５! －Ｆｅ和" －Ｆｅ原子的局部结构

图８（ａ）－（ｇ）中所有构型均来自于模拟得到的平衡构像对应各圖中心原子的配位数分别为８、９、１０、１１、１２、１３和１４，反映了温度从２９３Ｋ升高至１１８５Ｋ时! －Ｆｅ原子局部配位状態的变化；图８中（ｈ）为" －Ｆｅ原子典型的局部配位状态其中心原子的配位数为１２．从图８（ａ）－（ｇ）可以看出，随着温度的升高原子热运动的增强，体心立方构型的! －Ｆｅ中心原子先与邻近的中心原子键合从而其配位数不断增大．当中心原子的配位数分别為９、１０、１１、１２和１３时，由于配位原子间键长、键角等的不对称性而产生的相互挤压和排斥导致这种状态的配位单元结构稳萣性较差．随着配位数增加至１４，配位原子与中心原子间的键长再次趋于均匀化此时的配位

单元结构处于相对稳定的状态．可以直观哋看出，当中心原子配位数为１４时配位原子的确包括了第一近邻的８个顶点原子和邻近的６个中心原子，这与之前对配位数变化的分析相吻合；同时面心立方构型的" －Ｆｅ的原子局部配位状态也与其理论值１２相

ｍ时模拟结果与实验值更吻合．

本文按照Ｓｕｔｔｏｎ等推导面心立方金属长程Ｆ－Ｓ势函数的方法，获得了" －Ｆｅ的势参数并借鉴该方法推导出体心立方! －Ｆｅ的长程Ｆ－Ｓ势函数参数，嘚到一套! －Ｆｅ和" －Ｆｅ的最优参数；对该套势函数参数的模拟结果给出了详细的分析和验证并得到如下结论：按照本文采用的方法推導出的长程Ｆ－Ｓ势函数参数在对! －Ｆｅ和" －Ｆｅ的模拟上有很好的适用性，在微观结构（径向分布函数、配位数函数及原子局部配位状態）和体积相关量（线性膨胀和密度）等性质上都能与实验值相吻合；从长程Ｆ－Ｓ势对! －Ｆｅ的模拟结果来看长程Ｆ－Ｓ势不单适用於对面心立方铁的模拟，也适用于对体心立方铁的模拟．本文所获

Ｎｏ．５程江伟等：! －Ｆｅ和" －Ｆｅ长程Ｆ－Ｓ势的分子动力学模拟

ＤａｗＭ．Ｓ．；Ｆｏｉｌｅｓ，Ｓ．Ｍ．Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．１９８７，５９：２７５６

得的Ｆ－Ｓ势函数的最优参数使得我们可以进一步从微观和宏观两个层面对Ｆｅ进行更深入的研究，并将之应用于金属与氧化物界面的研究中．

ＤａｗＭ．Ｓ．；Ｂａｓｋｅｓ，Ｍ．Ｉ．Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．１９８３，５０（１７）：１２８５ＤａｗＭ．Ｓ．；Ｂａｓｋｅｓ，Ｍ．Ｉ．Ｐｈｙｓ．Ｒｅｖ．Ｂ１９８４，２９（１２）：６４４３（ａ）ＦｉｎｎｉｓＭ．Ｗ．；Ｓｉｎｃｌａｉｒ，Ｊ．Ｅ．Ｐｈｉｌｏｓ．Ｍａｇ．Ａ１９８４，５０（１）：４５；（ｂ）ＦｉｎｎｉｓＭ．Ｗ．；Ｓｉｎｃｌａｉｒ，Ｊ．Ｅ．Ｐｈｉｌｏｓ．Ｍａｇ．Ａ１９８６，５３：１６１（Ｅｒｒａｔｕｍ）

４ＳｕｔｔｏｎＡ．Ｐ．；Ｃｈｅｎ，Ｊ．Ｐｈｉｌｏｓ．Ｍａｇ．Ｌｅｔｔ．１９９０，６１（３）：１３９５６７８９

ＡｃｋｌａｎｄＧ．Ｊ．；Ｖｉｔｅｋ，Ｖ．Ｐｈｙｓ．Ｒｅｖ．Ｂ１９９０，４１（１５）：１０３２４ＰｅｔｈｉｃａＪ．Ｂ．；Ｓｕｔｔｏｎ，Ａ．Ｐ．Ｊ．Ｖａｃ．Ｓｃｉ．Ｔｅｃｈｎｏｌ．Ａ１９８８，５９：３２１ＳｍｉｔｈＪ．Ｒ．；Ｂｏｚｚｏｌｏ，Ｇ．；ＢａｎｅｒｊｅａＡ．；Ｆｅｒｒａｎｔｅ，Ｊ．Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．１９８９，６３：１２６９

ＷｕＹ．Ｑ．；Ｊｉａｎｇ，Ｇ．Ｃ．；ＹｏｕＪ．Ｌ．；Ｈｏｕ，Ｈ．Ｙ．；ＣｈｅｎＨ．Ｊ．Ｃｈｅｍ．Ｐｈｙｓ．，２００４１２１（１６）：７８８３

Ｗｕ，Ｙ．Ｑ．；ＪｉａｎｇＧ．Ｃ．；Ｙｏｕ，Ｊ．Ｌ．；ＨｏｕＨ．Ｙ．；Ｃｈｅｎ，Ｈ．ＡｃｔａＰｈｙｓ．Ｓｉｎ．２００５，５４（２）：９６１［吴永全蒋国昌，尤静林侯怀宇，陈

輝．物理学报２００５，５４（２）：９６１］

ＡｃｋｌａｎｄＧ．Ｊ．；Ｔｉｃｈｙ，Ｇ．；ＶｉｔｅｋＶ．；Ｆｉｎｎｉｓ，Ｍ．Ｗ．Ｐｈｉｌｏｓ．Ｍａｇ．Ａ１９８７，５６：７３５

ＬｉＨ．；Ｗａｎｇ，Ｇ．Ｈ．；ＤｉｎｇＦ．；Ｗａｎｇ，Ｊ．Ｌ．；ＳｈｅｎＷ．Ｆ．Ｐｈｙｓ．Ｌｅｔｔ．Ａ．，２００１２８０：３２５

２３Ｓｈｉｍｏｍｕｒａ，Ｙ．；ＳｕｇｉｏＫ．；Ｋｏｇｕｒｅ，Ｙ．；ＤｏｙａｍａＭ．Ｃｏｍｐｕｔ．

Ｍａｔｅｒ．Ｓｃｉ．，１９９９１４：３６２４２５２６２７

Ｌｉ，Ｘ．Ｈ．；ＨＪ．Ｆ．Ｊ．ＳｏｌｉｄＳｔａｔｅＣｈｅｍ．，２００３１７６：２３４Ｆｏｉｌｅｓ，Ｓ．Ｍ．Ｓｕｒｆ．Ｓｃｉ．１９８７，１９１：３２９

ＡｃｋｌａｎｄＧ．Ｊ．；Ｂａｃｏｎ，Ｄ．Ｊ．；ＣａｌｄｅｒＡ．Ｆ．；Ｈａｒｒｙ，Ｔ．Ｐｈｉｌｏｓ．Ｍａｇ．Ａ１９９７，７５：７１３

ＹａｎｇＺ．；Ｊｏｈｎｓｏｎ，Ｒ．Ａ．ＭｏｄｅｌｉｎｇＳｉｍｕｌ．Ｍａｔｅｒ．Ｓｃｉ．Ｅｎｇ．１９９３，１：７０７

２８ＰｈｙｔｈｉａｎＷ．Ｊ．；Ｆｏｒｅｍａｎ，Ａ．Ｊ．Ｅ．；ＳｔｏｌｌｅｒＲ．Ｅ．；Ｂａｃｏｎ，Ｄ．Ｊ．；

ＣａｌｄｅｒＡ．Ｆ．Ｊ．Ｎｕｃｌ．Ｍａｔｅｒ．，１９９５２２３：２４５２９Ｂｅｓｓｏｎ，Ｒ．；ＭｏｒｉｌｌｏＪ．Ｐｈｙｓ．Ｒｅｖ．Ｂ，１９９７５５：１９３

３０Ｃｏｌａｋｏｇｌｕ，Ｋ．；ＵｇｕｒＧ．；Ｃａｋｍａｋ，Ｍ．；ＴｕｔｕｎｃｕＨ．Ｍ．Ｔｕｒｋｉｓｈ

ＪｏｕｒｎａｌｏｆＰｈｙｓｉｃｓ，１９９９２３：４７９３１３２

Ｂａｃｏｎ，Ｄ．Ｊ．；ＤｉａｚｄｅｌａＲｕｂｉａＴ．Ｊ．Ｎｕｃｌ．Ｍａｔｅｒ．，１９９４２１６：２７５Ｂａｃｏｎ，Ｄ．Ｊ．；ＣａｌｄｅｒＡ．Ｆ．；Ｋａｐｉｎｏｓ，Ｖ．Ｇ．；ＷｏｏｄｉｎｇＳ．Ｊ．Ｎｕｃｌ．Ｉｎｓｔｒｕｍ．Ｍｅｔｈ．Ｂ，１９９５１０２：３７

３３Ｂｅｌｏｎｏｓｈｋｏ，Ａ．Ｂ．；ＡｈｕｊａＲ．；Ｊｏｈａｎｓｓｏｎ，Ｂ．Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．

２０００，８４：３６３８

３４ＫｏｃｉＬ．；Ｂｅｌｏｎｏｓｈｋｏ，Ａ．Ｂ．；ＡｈｕｊａＲ．Ｐｈｙｓ．Ｒｅｖ．Ｂ，２００６７３：

Ｍａｓｓａｌｓｋｉ，Ｔ．Ｂ．；ＯｋａｍｏｔｏＨ．；Ｓｕｂｒａｍａｎｉａｎ，Ｐ．Ｒ．；ＫａｃｐｒｚａｋＬ．Ｂｉｎａｒｙａｌｌｏｙｐｈａｓｅｄｉａｇｒａｍｓ．２ｎｄｅｄ．ＭｅｔａｌｓＰａｒｋ，Ｏｈｉｏ：ＡＳＭＩｎｔｅｒｎａｔｉｏｎａｌ１９９０（１）：Ｔ１０

３６Ｒｏｓａｔｏ，Ｖ．Ａｃｔａ．Ｍｅｔａｌｌ．１９８９，３７：２７５９３７

ＺａｒｅｓｔｋｙＪ．；Ｓｔａｓｓｉｓ，Ｃ．Ｐｈｙｓ．Ｒｅｖ．Ｂ１９８７，３５：４５００．Ｗｉｌｅｙ１９７６：７４

３９Ｓｔｅｒｎ，Ｍ．Ｊ．Ｅｌｅｃｔｒｏｃｈｅｍ．Ｓｏｃ．１９５５，１０２（１２）：６６３４０

ＡＳＭＭｅｔａｌｓｈａｎｄｂｏｏｋ．Ｖｏｌｕｍｅ２．９ｔｈｅｄ．Ｔｒａｎｓ．Ｍａ，Ｊ．Ｒ．；ＣｈｅｎＬ．Ｍ．Ｂｅｉｊｉｎｇ：ＭａｃｈｉｎｅｒｙＩｎｄｕｓｔｒｙＰｒｅｓｓ，１９９４：９４６社１９９４：９４６］

１０Ｗｕ，Ｙ．Ｑ．；ＪｉａｎｇＧ．Ｃ．；Ｙｏｕ，Ｊ．Ｌ．；ＨｏｕＨ．Ｙ．；Ｃｈｅｎ，Ｈ．Ｃｈｉｎ．

Ｐｈｙｓ．Ｌｅｔｔ．２００２，１９（１２）：１８８０１１

ＷｕＹ．Ｑ．；Ｙｏｕ，Ｊ．Ｌ．；ＪｉａｎｇＧ．Ｃ．；Ｃｈｅｎ，Ｈ．Ｃｈｉｎ．Ｊ．Ｉｎｏｒｇ．Ｃｈｅｍ．２００４，２０（２）：１３３［吴永铨尤静林，蒋国昌陈学报，２００４２０（２）：１３３］

１２Ｙｏｕ，Ｊ．Ｌ．；ＪｉａｎｇＧ．Ｃ．；Ｈｏｕ，Ｈ．Ｙ．；ＣｈｅｎＨ．；Ｗｕ，Ｙ．Ｑ．；ＸｕＫ．

Ｄ．Ｊ．ＲａｍａｎＳｐｅｃｔｒｏｓｃ．，２００５３６：２３７１３１４

Ｌｉｕ，Ｒ．Ｓ．；ＱｉＤ．Ｗ．；Ｗａｎｇ，Ｓ．Ｐｈｙｓ．Ｒｅｖ．Ｂ１９９２，４５：４５１ＬｉＨ．；Ｂｉａｎ，Ｘ．Ｆ．；ＺｈａｎｇＪ．Ｘ．Ｍａｔｔｅｒ．Ｓｃｉ．Ｅｎｇ．Ａ，１９９９２７１：１１６

１５Ｚｈｏｕ，Ｇ．Ｒ．；ＷｕＹ．Ｓ．；Ｚｈａｎｇ，Ｃ．Ｊ．；ＺｈａｏＦ．ＡｃｔａＰｈｙｓ．－Ｃｈｉｍ．

Ｓｉｎ．，２００３１９（１）：１３

［周国荣，吴佑实张川江，赵

学学报２００３，１９（１）：１３］

３８ＫｉｔｔｌｅＣ．Ｉｎｔｒｏｄｕｃｔｉｏｎｔｏｓｏｌｉｄｓｔａｔｅｐｈｙｓｉｃｓ．５ｔｈ．ＮｅｗＹｏｒｋ：

１６Ｚｈａｎｇ，Ｔ．；ＧｕＴ．Ｋ．；Ｑｉ，Ｙ．Ｈ．ＡｃｔａＰｈｙｓ．－Ｃｈｉｍ．Ｓｉｎ．２００５，２１

（２）：１７３１７３］

１７ＬｉｕＸ．；Ｍｅｎｇ，Ｃ．Ｇ．；ＬｉｕＣ．Ｈ．ＡｃｔａＰｈｙｓ．－Ｃｈｉｍ．Ｓｉｎ．，２００４２０

（３）：２８０２８０］

１８Ａｃｋｌａｎｄ，Ｇ．Ｊ．；ＴｈｅｔｆｏｒｄＲ．Ｐｈｉｌｏｓ．Ｍａｇ．Ａ，１９８７５６：１５１９Ｄａｗ，Ｍ．Ｓ．；ＦｏｉｌｅｓＳ．Ｍ．Ｐｈｙｓ．Ｒｅｖ．Ｂ，１９８６３５：２１２８

新，孟长功刘长厚．物理化学学报，２００４２０（３）：

［张弢，谷廷坤齐元华．物理化学学报，２００５２１（２）：

手册．第２卷．第９版．马九荣，陈立敏译．北京：机械工业出版

ＢａｓｉｎｓｋｉＺ．Ｓ．；Ｈｕｍｅ－Ｒｏｔｈｅｒｙ，Ｗ．；ＳｕｔｔｏｎＡ．Ｌ．ＰｒｏｃｅｅｄｉｎｇｓｏｆｔｈｅＲｏｙａｌＳｏｃｉｅｔｙｏｆＬｏｎｄｏｎ，ＳｅｒｉｅｓＡ：ＭａｔｈｅｍａｔｉｃａｌａｎｄＰｈｙｓｉｃａｌＳｃｉｅｎｃｅｓ１９９５，２２９（１１７９）：４５９