石墨烯力学性能的分子动力学模拟：周期性边界与自由边界对断裂行为的系统性影响

石墨烯的力学传奇从2008年Lee等人用纳米压痕实验测出约130GPa的断裂强度开始，分子动力学模拟界就陷入了一场持续至今的方法论拉锯。周期性边界条件（PBC）下的模拟结果总是漂亮得令人不安——单轴拉伸沿armchair方向，应力-应变曲线一路攀升到110GPa才出现那道标志性的断裂陡坡，对应的断裂应变稳稳停在22%。换成自由边界条件（FBC）重新跑一遍，峰值应力直接跌到78GPa，断裂应变缩水到14%。30%的强度差异不是数值噪声，而是边界条件在微观尺度上重新分配应力后的必然结果。

这种差异的根源藏在边缘那两排碳原子里。PBC下，石墨烯薄片在x和y方向被无限复制，每条模拟盒子的”边界”实际上都是体相区域，应力可以在虚拟的周期性映像之间连续传递。FBC则把边缘原子赤裸裸地暴露出来——悬键碳原子缺少相邻晶格的约束，局部配位数从理想的3降到了2甚至1，这些边缘原子在应变加载初期就承受了远超体相原子的局部应力。OVITO里把局部von Mises应力映射成彩色云图，边缘两排原子总是最先被染成深红色，像刀切一样整齐地沿着裂纹路径排列。

AIREBO与ReaxFF对C-C键断裂的描述分歧

势函数的选择在这里不是一个小问题。AIREBO势（Adaptive Intermolecular Reactive Empirical Bond-Order）对石墨烯C-C键的键级描述基于Tersoff-Brenner框架，键的断裂由一个平滑的键级函数控制，在键长超过约1.85Å时键级快速衰减到零。这种处理方式在石墨烯拉伸模拟中给出的断裂行为干净利落——应力达到峰值后，边缘某处的一个C-C键突然断裂，应力集中瞬间向两侧传播，整个薄片在皮秒量级内完成脆断。

ReaxFF则走了一条完全不同的路。作为反应力场，它在每一步动力学积分中都会重新计算键级，电荷分布也会随原子位置动态更新。这种”实时感知”让ReaxFF在描述石墨烯断裂时呈现出更丰富的细节——边缘悬键处的电荷重新分布会轻微改变局部键能，使得断裂路径不再是一条直线，而是沿着zigzag或armchair边缘的晶向蜿蜒前进。代价是计算量，同样规模的石墨烯薄片（约5000个原子），ReaxFF的每一步积分耗时是AIREBO的8-12倍。在同等计算预算下，AIREBO可以跑到20ns的拉伸过程，ReaxFF往往只能覆盖前5ns。

两种势函数在PBC下的断裂强度预测倒是没有本质分歧——AIREBO给出108-112GPa，ReaxFF给出103-109GPa，差异在5%以内。真正的分歧出现在FBC下对断裂路径的预测。AIREBO倾向于让裂纹沿初始缺陷垂直传播，ReaxFF则会因为电荷介导的局部键能波动而产生分叉裂纹，后者与实验观测到的石墨烯断裂形貌吻合得更好。

孔洞缺陷对边界效应的放大机制

干净的石墨烯薄片在PBC和FBC之间的差异固然显著，引入孔洞缺陷后，这种差异会被进一步放大到一个无法忽视的程度。在薄片中心放置一个直径1nm的圆形孔洞（约移除30个碳原子），PBC下的断裂强度从110GPa降到89GPa，降幅约19%。FBC下则从78GPa跌到54GPa，降幅达到31%。孔洞像一面放大镜，把边界条件带来的应力分配差异投射到了整个薄片的应力场中。

孔洞直径从1nm增加到3nm的过程中，两种边界条件的断裂强度差距从34GPa扩大到42GPa。3nm孔洞在石墨烯晶格中移除了约270个碳原子，孔洞边缘的应力集中因子在弹性理论中可以估算为3——这意味着远场应力为30GPa时，孔洞边缘的局部应力已经突破90GPa。FBC下的边缘悬键区域和孔洞边缘的应力集中区如果距离较近（小于5nm），两者会产生叠加效应，局部应力峰值可以超出远场应力4-5倍，这是PBC下永远不会出现的场景，因为PBC没有”边缘悬键区域”这个概念。

应变方向的各向异性响应

石墨烯的六方晶格赋予了它强烈的各向异性，armchair方向和zigzag方向在力学响应上的差异在边界条件切换时表现出令人玩味的不对称性。沿armchair方向加载单轴应变，PBC下断裂强度110GPa，FBC下78GPa，差距32GPa。换成zigzag方向，PBC下断裂强度升至约121GPa（源于C-C键与加载方向的几何关系），FBC下降到86GPa，差距35GPa。armchair方向对边界条件的变更更敏感——这不是因为armchair方向的本征强度更低，而是因为armchair晶向在FBC下的边缘原子配位数下降更严重，悬键效应更突出。

这种各向异性在LAMMPS输入脚本中的体现非常直接——fix deform命令的dir参数设为x还是y，对应的晶向定义取决于初始建模时如何用lattice custom命令定义石墨烯的原胞。很多模拟研究中断裂强度的”异常分散”实际上源于建模时晶向定义的模糊，而非势函数或边界条件本身的问题。建议在模拟开始前用write_data输出初始结构，在OVITO中确认晶向与预期一致，这个小步骤可以省去后续大量关于”为什么我的断裂强度只有70GPa”的排查时间。

模拟参数对结果的隐性控制

除了边界条件和势函数，若干看似次要的模拟参数实际上对断裂行为有系统性影响。拉伸速率是最容易被低估的变量——文献中常见的应变率范围是0.001/ps到0.01/ps，在这个区间内，断裂强度随应变率的对数近似线性增长，斜率约为每数量级8-12GPa。将应变率从0.001/ps降到0.0001/ps（需要更长的模拟时间），断裂强度会进一步下降5-8GPa，这个趋势与实验观测一致：更慢的加载速率允许更多的热涨落参与键断裂的成核过程，使得断裂在更低的远场应力下发生。

温度控制方案也有发言权。fix npt与fix temp/rescale给出的断裂强度可以相差15%以上，前者允许应力在三个方向上自洽调整，后者强制维持温度可能导致应力场的非物理涨落。石墨烯薄片的面内热膨胀系数在300K附近接近于零甚至为负，这意味着在NPT系综下，面内应力可能自发产生约0.1GPa的热应力，对于追求高精度断裂强度的工作，这个量级已经不可忽略。

结论与模拟策略建议

边界条件的选择不应是一个默认选项，而是需要与研究目标主动对齐的决策。研究石墨烯本征力学性能（如理论强度上限）应优先使用PBC，并将结果与无孔、无缺陷的理想晶格对应。研究石墨烯纳米带的边缘效应、裂纹传播路径或与实际器件几何对齐时，FBC是不可回避的选择，此时建议在薄片宽度上取三个以上不同值做收敛性测试——当宽度超过10nm时，FBC下的断裂强度会逐渐向PBC结果靠拢，因为边缘效应在绝对尺度上被稀释了。

势函数的选择上，AIREBO适合大尺度（>10000原子）、长时间（>10ns）的筛选性模拟；ReaxFF适合小尺度但需要精确描述键断裂化学细节的场景。孔洞缺陷的引入会放大两种边界条件的差异，在含缺陷样品的模拟中，建议同时报告PBC和FBC下的结果，并明确标注边缘到缺陷的距离，为读者提供完整的应力场上下文。

欢迎访问 keyanxueshu.com 了解分子动力学计算服务。

参考文献

1. Lee, C., Wei, X., Kysar, J. W., & Hone, J. (2008). Measurement of the elastic properties and intrinsic strength of monolayer graphene. Science, 321(5887), 385-388.
2. Stuart, S. J., Tutein, A. B., & Harrison, J. A. (2000). A reactive potential for hydrocarbons with intermolecular interactions. Journal of Chemical Physics, 112(14), 6472-6486.
3. Liu, P., Zhang, Y., & Cheng, Y. (2011). Detection of initial damage in graphene nanoribbons under uniaxial tension by MD simulations. Computational Materials Science, 50(8), 2515-2521.
4. Brenner, D. W., et al. (2002). A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. Journal of Physics: Condensed Matter, 14(4), 783-802.