分子动力学模拟揭秘：非晶材料断裂韧性的原子尺度起源

1. 项目概述：从原子视角窥探玻璃的“筋骨”

玻璃，这种我们日常生活中再熟悉不过的材料，其本质是一种非晶态固体。与金属或晶体不同，玻璃内部的原子排列没有长程有序的周期性，呈现出一种“冻结的液体”状态。这种无序性赋予了玻璃透明、各向同性等独特性质，但也使其在力学性能上，尤其是断裂行为，充满了复杂性和不确定性。为什么有些玻璃一碰就碎，而有些却能承受巨大的弯曲？这个问题的答案，深埋在原子尺度的微观世界里。

传统的材料研究多依赖于宏观实验，但面对非晶材料这种“黑箱”，宏观数据往往难以揭示其内在的物理机制。分子动力学模拟，就像一台超级显微镜，允许我们追踪每一个原子的运动轨迹，在计算机中“制备”出具有特定微观结构的材料模型，并对其施加虚拟的载荷，从而直观地观察从弹性变形到最终断裂的全过程。本次分享的项目，正是运用这一利器，聚焦于一个经典的模型体系——二维二氧化硅玻璃。选择二维体系并非为了简化，而是为了剥离三维空间的复杂性，更清晰地聚焦于非晶结构本身对力学行为的影响。我们构建了三种具有代表性的微观结构模型：纳米晶玻璃、准晶玻璃和连续随机网络玻璃，旨在系统性地回答一个核心问题：原子尺度上的无序结构差异，如何一步步决定材料的宏观断裂韧性？

对于材料科学、物理或相关领域的从业者，无论是初涉模拟的研究生，还是希望深化对非晶材料理解的技术人员，这篇文章将带你深入一场从原子建模到宏观性能解读的完整旅程。我们将不仅复现模拟的关键步骤，更会拆解每一个结果背后的物理图像，分享参数设置背后的考量，以及如何从海量的原子轨迹数据中提炼出有价值的科学见解。

2. 核心模型构建与结构表征解析

分子动力学模拟的第一步，也是决定模拟成败与物理意义的关键一步，就是构建能够真实反映材料特性的原子模型。对于二维二氧化硅，其基本结构单元是硅氧四面体在二维平面上的投影，即一个硅原子与三个氧原子键合形成的三角形单元，这些单元通过共用氧原子连接成网络。

2.1 三种非晶结构模型的生成策略

我们构建了三种截然不同的二维二氧化硅玻璃模型，它们代表了非晶结构中从“有序”到“无序”的谱系。

2.1.1 纳米晶玻璃模型顾名思义，NCG模型由许多微小的晶粒镶嵌在非晶基质中构成。我们的生成方法是：首先，在模拟盒子中随机撒播多个晶核点，每个晶核以二维β-鳞石英的晶体结构（一种二氧化硅的晶体相）为模板进行生长。通过控制晶核的数量和生长速率，我们可以获得不同晶粒尺寸和体积分数的纳米晶结构。当晶粒生长到相互接触或达到预设尺寸后，剩余的空间则通过快速淬火（即高温熔体快速冷却）的方式形成非晶基质。这个过程模拟了某些通过控制结晶工艺制备的玻璃陶瓷材料。

注意：晶粒与基质的界面是关键。在模拟中，需要确保界面处的原子通过能量最小化进行了充分弛豫，以避免引入不切实际的高能界面，从而影响断裂起始位置的判断。

2.1.2 准晶玻璃模型Para模型比NCG更加无序，但保留了中程有序的特征。我们采用了一种“拓扑约束”生长算法。首先生成一个完全随机的初始原子构型，然后通过引入一组约束条件（如键长、键角的目标分布函数，这些函数来源于对实验或高质量CRN结构的统计）进行迭代优化。在优化过程中，系统被允许形成一些局域的、类晶的短程有序团簇，但这些团簇之间没有明确的晶界，呈现出一种“近程有序、中程弥散”的状态。这类似于某些通过特殊工艺制备的、结构高度均匀的非晶合金。

2.1.3 连续随机网络模型CRN模型是描述理想玻璃态最经典的模型。它假设网络是连续、无限且完全无序的，但每个原子都满足确定的配位数（硅为4，氧为2）。我们采用Wooten-Winer-Weaire算法或其变体来生成。该算法从一个晶态种子开始，通过一系列随机的键切换操作来引入无序，同时严格保持网络的连通性和配位数约束，直到系统的能量和径向分布函数收敛到典型玻璃态的特征。最终得到的结构没有任何周期性，是三种模型中最无序的。

2.2 关键结构表征量的计算与物理意义

模型建好后，如何定量描述它们的结构差异？这就需要借助一系列结构表征函数。

2.2.1 径向分布函数：探测短程有序的“尺子”RDF是材料科学中最基础的结构表征工具。对于我们的体系，我们计算了Si-Si、Si-O和O-O的RDF。计算原理是：以某个原子为中心，统计在距离r到r+dr的球壳内发现其他原子的平均数量，然后对体系中的所有原子和所有构型进行平均并归一化。

• 第一峰位置：对应最近邻原子间的平均距离。例如，Si-O RDF的第一峰位置直接给出了Si-O键的平均键长，这是判断势函数合理性和结构弛豫程度的重要指标。
• 第一峰面积：积分后可以得到配位数。Si-O第一峰面积应接近4（Si的配位数），O-Si第一峰面积应接近2（O的配位数）。
• 峰的宽度和衰减：峰越尖锐，说明键长分布越集中，有序度越高。从第二、第三峰开始，对于晶体，RDF会出现清晰的峰；对于CRN玻璃，这些峰会迅速展宽并湮灭，这正是长程无序的体现。将模拟CRN的RDF与实验测得的RDF（如图1e,f）进行对比，是验证模型可靠性的黄金标准。

2.2.2 环分布统计：揭示中程有序的“拓扑地图”在二维网络中，环（由Si和O原子交替连接形成的闭合多边形）是描述拓扑结构的核心。我们统计了网络中所有最小环的边数（通常是4-10元环）。对于完美的β-鳞石英晶体，只存在6元环。而在非晶网络中，会出现4、5、7、8元环等多种拓扑缺陷。

• 环分布图：如图1d所示，NCG模型由于含有晶粒，其环分布会在6元环处出现一个尖锐的峰，同时在其他位置有非晶基质贡献的宽分布。Para模型的分布较宽，但可能在6元环附近仍有凸起。CRN模型的分布最宽且平缓，是典型的无特征分布。环分布的差异直接影响了网络的刚性、自由体积分布，进而影响力学性能。

2.2.3 角度分布函数与结构相似性参数ADF统计了诸如Si-O-Si和O-Si-O等键角的分布。对于二氧化硅，Si-O-Si键角分布通常在120°-180°之间，中心约在145°。ADF的宽度反映了网络键角扭曲的程度。结构相似性参数s则是一个更全局的量度，用于量化局部原子环境与某个参考结构（如理想的晶体环境）的差异，其概率分布函数（图1g）可以直观显示体系内结构异质性的程度。

通过这些细致的结构表征，我们不仅在视觉上（图1a-c），更在定量上严格区分了三种模型，为后续关联结构与性能奠定了坚实的基础。模拟得到的CRN的RDF和ADF与实验数据高度吻合（图1e,f,h,i），这让我们对模拟的可靠性充满信心。

3. 分子动力学模拟流程与关键参数设置

有了可靠的原子模型，下一步就是让它们“动”起来，接受力学测试。分子动力学模拟的本质是数值求解牛顿运动方程，其流程和参数设置直接决定了模拟的物理真实性和计算效率。

3.1 势函数选择：描述原子间相互作用的“宪法”

对于二氧化硅体系，原子间的相互作用通常由经验势函数来描述。我们选择了基于BKS模型的修正势函数。BKS势是一种经典的二体势，包含了长程库仑相互作用和短程Buckingham势的形式。

E_ij = q_i q_j / r_ij + A_ij exp(-B_ij r_ij) - C_ij / r_ij^6

其中，第一项是库仑项（q_i为原子电荷），第二项和第三项是短程排斥和吸引项。我们采用Wolf方法对长程库仑力进行截断和修正，以在保证精度的前提下大幅提升计算速度。

实操心得：势函数的选择是分子动力学模拟的基石。对于二氧化硅，除了BKS，还有TTAM、CHIK等势函数。选择时需权衡精度与速度。BKS势在描述玻璃态结构（如RDF、密度）方面表现良好，且计算效率较高，适合进行大规模断裂模拟。务必在正式模拟前，用小体系测试所选势函数能否复现实验已知的键长、键角、弹性模量等基本性质。

3.2 系统弛豫与平衡态获取

在施加载荷前，必须让初始构建的模型达到平衡态，即体系能量和压力等宏观量在平均值附近小幅波动。

1. 能量最小化：首先对初始模型进行共轭梯度法能量最小化，消除原子间的过度重叠和不合理的高应力。
2. NPT系综弛豫：在恒温恒压系综下运行一段时间（通常数十皮秒），让体系密度弛豫到平衡值。温度通过Nosé-Hoover热浴控制，压力通过Parrinello-Rahman控压方法控制。对于二维体系，需注意控制面内压力，而垂直于平面方向（第三维）的压力设为0。
3. NVT系综弛豫：在平衡体积下，切换至恒温恒容系综继续弛豫，确保温度均匀，体系处于稳定的玻璃态。我们通常将温度设定在远低于玻璃化转变温度Tg的数值（如300K）。

关键参数：弛豫时间需要足够长。我们通过监测体系总能量、温度、压力以及RDF是否不再有漂移来判断平衡。对于含有数千原子的体系，NPT+NVT弛豫总计需要至少100-200皮秒。

3.3 单轴拉伸模拟的实现细节

我们采用位移控制的方式进行单轴拉伸模拟，这是研究断裂行为的常用方法。

1. 加载方向：沿模拟盒子的x方向（armchair方向）施加应变。同时，允许y方向自由收缩（泊松效应），而z方向（二维平面外）保持固定厚度。
2. 应变率：这是一个至关重要的参数。分子动力学模拟受限于计算资源，使用的应变率（通常为10^8 - 10^9 /s）远高于实验速率（10^-3 - 10^0 /s）。如此高的应变率会抑制原子的热激活运动，可能导致材料表现得比实际更脆。为了部分抵消这一效应，我们采用了相对“慢”的应变率（在MD尺度下），例如5×10^8 /s，并通过与不同应变率的对比测试，评估应变率效应的影响。
3. 模拟步骤：在每个模拟时间步（通常为1飞秒），将模拟盒子在x方向的长度增加一个微小量ΔL = ε˙ * L_x0 * Δt，其中ε˙是应变率，L_x0是初始长度，Δt是时间步长。然后重新标定所有原子的x坐标（仿射变形），接着在NVT系综下运行若干步让原子弛豫到新的构型，并计算维里应力。
4. 应力计算：体系在变形过程中的应力张量通过维里定理计算，包含了动能项和原子间相互作用力的贡献。我们主要关注沿拉伸方向的轴向应力（σ_xx）。

关键技巧：为了捕捉断裂起始的精确位置，我们在整个拉伸过程中，以很高的频率（如每100步）输出原子的位置、速度和应力数据。这些海量的轨迹文件是后续所有分析的源头。

4. 断裂起始机制：从最弱环节到空洞形核

材料不会无缘无故断裂，断裂总是从最薄弱的环节开始。通过分子动力学模拟，我们可以像侦探一样，在断裂发生前就锁定这些“嫌疑位置”，并观察“犯罪现场”是如何形成的。

4.1 维里系数：预测断裂起始点的微观应力探针

宏观应力是体系整体的平均值，而断裂往往由局部极高的应力集中引发。维里系数V_θ是我们定义的一个局部应力度量。对于每个原子i，我们计算其在一个局部区域（例如其最近邻壳层）内的应力张量，然后通过坐标变换，得到沿拉伸方向（θ）的正应力分量V_θ(i)。这个值直观反映了该原子所在局部区域被拉伸的紧张程度。

如图2a所示，在未变形的平衡结构中，我们就可以绘制出V_θ的空间分布图。那些呈现高正值（红色区域）的原子团簇，意味着即使在零应变下，该处的局部网络也处于一种“预拉伸”的亚稳态。它们就是潜在的断裂起始点。我们的模拟结果惊人地显示，最终发生空洞形核的位置（图2a中圆圈标注），与这些高V_θ区域高度重合。

4.2 空洞形核的临界条件与结构关联

我们统计了所有模拟中空洞形核时的应变值（ε_nuc）。分析发现，空洞形核应变与平衡结构中的两个关键结构参数存在强烈的相关性：

1. 与最拉伸状态的关系：对于每个样本，我们找到其平衡结构中V_θ的最小值（即最拉伸的状态）。如图2b所示，V_θ_min与ε_nuc之间存在近似指数衰减的关系。V_θ_min越小（即初始就越拉伸），材料越早发生空洞形核，韧性越差。这为预测材料的断裂起始韧性提供了一个基于初始结构的微观判据。
2. 与Si-O-Si键长的关系：我们测量了所有Si-O-Si桥氧键的长度（即两个硅原子之间的距离）。统计发现，空洞往往在那些具有较长Si-O-Si键的附近形核（图2d）。较长的键通常意味着键能较弱，或者该处的网络拓扑结构较为松散（自由体积较大），在拉伸下更容易被拉开。

注意事项：空洞形核是一个统计过程。即使在同一模型中，由于原子级结构的随机涨落，每次模拟的形核位置和应变也可能不同。因此，我们需要对每种结构进行大量（如50-100次）的独立模拟，并统计其累积概率分布函数（CDF），如图2c所示。CRN的CDF最宽，说明其形核应变的波动性最大，这与它结构的均匀性（没有明显的优先形核点）有关；而NCG的CDF较窄，因为晶界作为明确的缺陷，总是优先形核。

4.3 断裂事件的原子尺度快照

通过分析断裂瞬间的原子构型（图3d, g, j），我们可以清晰地看到空洞形核的微观过程：

• 在NCG中，空洞几乎无一例外地在晶界处形核。晶界处原子排列不规则，键合较弱，是天然的弱点。
• 在Para中，空洞在那些中程有序被破坏的局部松散区域形核。
• 在CRN中，空洞形核位置看起来更随机，但追溯其平衡结构，总能找到对应的V_θ高点或长键位置。

形核后，空洞并非立即导致灾难性断裂。它首先会稳定生长，直到达到一个临界尺寸。

5. 裂纹扩展动力学与韧性差异分析

空洞形核只是断裂的序幕，随后的裂纹扩展阶段才真正决定了材料的断裂韧性和断裂形貌。这是三种结构模型表现出最大差异的环节。

5.1 应力-应变响应与宏观韧性指标

图3a展示了三种玻璃和完美晶体在单轴拉伸下的应力-应变曲线。这是最宏观的性能对比：

• 晶体：表现出典型的脆性断裂特征——线弹性变形直至应力峰值，随后应力骤降，断裂应变小。
• NCG玻璃：曲线与晶体类似，但峰值应力（强度）和断裂应变略有降低。断裂是突然的，一旦裂纹从晶界萌生，便迅速扩展贯穿样品。
• Para玻璃：曲线出现了一定的非线性（塑性）特征，断裂应变明显增加。应力在达到峰值后并非瞬间跌落，而是有一个缓慢下降的过程。
• CRN玻璃：表现出最高的断裂应变和最显著的“延性”。应力-应变曲线在峰值后有一个很长的下降尾巴，意味着裂纹扩展过程是缓慢、受阻的。

峰值应力代表材料的强度，而应力-应变曲线下的面积则代表了材料断裂前吸收的能量，是衡量断裂韧性的直观指标。显然，CRN > Para > NCG ≈ Crystal。

5.2 裂纹扩展的两种模式与Griffith理论修正

我们追踪了裂纹尺寸（归一化为样品尺寸）随应变的变化（图4a,b）。发现裂纹扩展存在两种模式：

1. 直接快速扩展：主要发生在NCG和晶体中。一旦主裂纹形成，便迅速扩展，几乎不遇到阻力。
2. 受阻缓慢扩展伴随空洞化：主要发生在CRN和部分Para结构中。主裂纹前端会引发大量次级空洞的形核（图4c中void数量增加）。这些空洞或与主裂纹汇合（coalescence），或通过剪切带与主裂纹桥接（bridging），如图4e所示。这个过程消耗了大量额外能量，显著提高了韧性。

经典的Griffith理论认为，当裂纹尺寸达到临界值c_cr = 4γ / (πEε^2)时，裂纹将失稳扩展。其中γ是表面能，E是杨氏模量，ε是外加应变。我们从模拟中提取了γ和E，计算了c_cr（图4d中虚线）。结果显示，对于NCG，观测到的扩展裂纹尺寸与c_cr吻合较好。但对于CRN，许多裂纹在远小于c_cr时就开始缓慢扩展，而大量形核的次级空洞尺寸则小于c_cr且不扩展（非传播空洞）。这说明在非晶材料中，由于结构无序导致的应力场异质性（图4f），裂纹尖端并非处于理想的K场（应力强度因子场），局部波动很大。Griffith理论基于连续介质力学和均匀材料假设，在原子尺度的无序结构面前需要修正。

5.3 断裂路径粗糙度与“无序捕获”效应

我们定量分析了断裂后表面的粗糙度（图3b）。粗糙度定义为断裂轮廓线相对于平均高度的均方根值。

• NCG：断裂面最光滑，因为裂纹倾向于沿着相对平直的晶界扩展。
• CRN：断裂面最粗糙。裂纹在扩展过程中不断被局域的无序结构所偏折、分叉，甚至暂时“被困住”（disorder-trapping effect），需要更大的驱动力才能继续前进，如图4g所示。这种“之”字形的扩展路径极大地增加了新生表面积，从而消耗了更多能量（表面能γ本质上就是创造新表面所需的能量）。
• Para：粗糙度介于两者之间。

“空洞引导”的裂纹路径是CRN高韧性的另一个关键机制（图4g插图）。裂纹并非总是直线前进，而是倾向于朝着预先存在的或新形核的次级空洞方向扩展，因为连接两个空洞所需的能量可能低于直接撕裂完整的网络。这种路径选择进一步增加了裂纹扩展的曲折度。

6. 模拟实践中的常见问题与深度排查指南

分子动力学模拟非晶材料断裂是一个复杂的过程，从建模到分析，处处是坑。以下是我在大量实践中总结出的典型问题与解决方案。

6.1 模型构建与弛豫阶段问题

6.2 拉伸模拟与数据分析阶段问题

6.3 性能优化与计算资源管理

大规模分子动力学模拟极其消耗计算资源。一次包含数万原子、数百万时间步的模拟，在普通工作站上可能需要数周时间。

• 并行计算：务必使用支持并行（如MPI、OpenMP）的分子动力学软件，如LAMMPS。根据体系特点选择最优的并行分区策略。
• 邻居列表更新：合理设置邻居列表的截断半径和更新频率。对于断裂模拟，原子相对运动剧烈，需要比平衡模拟更频繁地更新邻居列表。
• 增量式应变加载：与其在每一步都重新标定所有原子坐标，不如采用“变形-弛豫”的循环。每次施加一个有限应变增量（如0.25%），然后在NVT下弛豫足够步数，再施加下一个增量。这更接近准静态过程，且允许在弛豫阶段使用更大的时间步长。
• 数据输出策略：全轨迹输出（所有原子的位置和速度）会产生TB级的数据。只输出你分析所必需的数据。例如，对于应力-应变曲线，只需输出系统整体的应力和应变；对于局部分析，可以只在高频输出特定原子组的信息，或者每隔很多步才输出一次全轨迹用于可视化。

通过这套系统的建模、模拟、分析和排错流程，我们才能从原子运动的噪声中，提炼出决定非晶材料断裂行为的清晰物理图像。这项研究不仅加深了对玻璃这一古老材料本质的理解，其揭示的“结构异质性调控韧性”的普适原理，也为设计更高性能的金属玻璃、高分子玻璃等非晶材料提供了宝贵的微观设计思路。