简介:本资源是面向材料科学与半导体物理研究者的开源动力学蒙特卡洛(KMC)模拟工具MMonCa的完整源代码包,专用于掺杂剂在晶体中的扩散行为建模及外延生长过程仿真,适用于研究生、科研人员及计算材料方向工程师开展机理研究与工艺参数优化。压缩包为ZIP格式,大小5.83MB,包含MMonCa-master主分支全部源码文件(含核心算法模块、示例输入配置、编译说明与基础文档),虽未提供具体文件总数与类型明细,但典型结构覆盖可编译C++源码、参数化输入模板及物理模型说明,便于本地部署与二次开发。已有568人学习下载,体现其在学术界的实际应用热度。用户可直接获取可运行的KMC框架,深入理解扩散跃迁概率设置、时间步长推进机制及表面吸附/脱附等关键物理过程的代码实现,同时基于开源特性灵活扩展新晶格模型或耦合热力学参数,显著提升材料生长模拟的自主建模能力。
1. MMonCa 是什么:一个专为晶体材料中掺杂原子扩散与外延生长建模而生的动力学蒙特卡洛(kMC)开源引擎
你有没有遇到过这样的困境:想模拟磷原子在硅单晶里怎么一步步跳到替位位置,或者砷掺杂剂在低温外延过程中如何被“冻”在间隙态、又在退火时突然活化——用第一性原理算单个跃迁势垒已经很吃力,而真实器件尺度下涉及百万级原子位点、跨越秒级时间尺度的演化,传统分子动力学(MD)早就算不动了。MMonCa 就是为解决这个「时间-空间尺度断层」而设计的:它不追踪每个原子的连续轨迹,而是把晶体抽象成一张带能量权重的跃迁图谱,用随机抽样+速率加权的方式,让掺杂原子在格点间“跳格子”,从而在合理计算开销下,复现扩散路径、偏析倾向、表面形核序列等关键动力学行为。它不是通用 kMC 框架(比如 Zacros 或 KMCLib),而是从晶体对称性、缺陷化学、表面台阶动力学出发深度定制的领域专用引擎。适合材料计算方向的研究生、半导体工艺仿真工程师,以及需要把 DFT 计算结果快速对接到器件级模型的团队——如果你手头已有掺杂原子的迁移能垒、结合能、表面吸附能等输入参数,MMonCa 就是你能把这些“纸面数据”变成“可演化的动态过程”的最短路径。
2. 从零跑通 MMonCa:编译、配置与第一个扩散模拟任务
MMonCa 的源码结构清晰,但它的构建逻辑和典型工作流与通用科学计算软件略有不同:它依赖晶体对称性解析和离散跃迁事件建模,因此必须先定义晶格拓扑,再注入物理参数,最后启动随机采样。整个流程不能跳步,否则连初始构型都生成不了。下面以硅(Si)中硼(B)掺杂的体扩散为例,带你走完最小可行闭环。
2.1 获取源码与环境准备:C++17 + CMake 3.16+ 是硬门槛
MMonCa 是纯 C++ 实现,无 Python 封装层,所有核心逻辑都在src/下。官方未提供预编译包,必须本地构建。我们推荐在 Linux 或 macOS 上操作(Windows 需 WSL2),确保系统已安装:
- GCC ≥ 10 或 Clang ≥ 12(需完整支持 C++17 并发与文件系统库)
- CMake ≥ 3.16(低于此版本无法识别
find_package(Threads)的现代写法) - Python 3.8+(仅用于辅助脚本,如
scripts/generate_lattice.py)
提示:不要用 conda 或 pip 安装的 CMake,它们常为旧版。直接从 cmake.org 下载二进制包并加入 PATH,或用
apt install cmake(Ubuntu 22.04+ 默认满足)。
克隆仓库后,进入根目录执行标准三步:
mkdir build && cd build cmake -DCMAKE_BUILD_TYPE=Release -DENABLE_MPI=OFF .. make -j$(nproc)-DENABLE_MPI=OFF是关键:MMonCa 当前不支持 MPI 并行,强行开启会导致链接失败。若你看到undefined reference to MPI_Init,一定是误启了 MPI 选项。编译成功后,build/bin/下会生成可执行文件mmonca。
2.2 构建硅晶格模板:用lattice_gen工具生成.lat文件
MMonCa 不接受任意坐标系输入,它要求晶格必须由lattice_gen工具生成标准.lat格式文件。该工具位于build/bin/lattice_gen,其输入是晶体学参数(非 POSCAR)。以金刚石结构硅为例(a = 5.431 Å,F d-3m 空间群),运行:
./lattice_gen \ --crystal=diamond \ --a=5.431 \ --output=si_4x4x4.lat \ --size=4,4,4 \ --vacancy=0.0 \ --interstitial=0.0参数说明:
--crystal=diamond:固定值,表示使用金刚石结构基元(Si、Ge 等均适用);--size=4,4,4:生成 4×4×4 个原胞的超胞,共 256 个格点(每个原胞含 8 个原子位点);--vacancy=0.0:初始空位浓度设为 0,避免干扰掺杂原子扩散主过程;--interstitial=0.0:同理,不预置间隙原子。
执行后生成si_4x4x4.lat,可用文本编辑器打开查看:前几行是晶格矢量与原胞内原子坐标(以分数坐标给出),随后是每个格点的类型标签(SIT表示替位位点,INT表示间隙位点)。这是 MMonCa 唯一认可的晶格描述方式——不能直接喂入 VASP 的 POSCAR 或 ASE 的 atoms 对象。
2.3 编写核心输入文件input.json:定义掺杂、跃迁与模拟策略
MMonCa 的行为完全由input.json控制。它不是自由格式 JSON,而是有严格 schema 的配置文件。以下是最小可用模板(对应 B 在 Si 中扩散):
{ "lattice_file": "si_4x4x4.lat", "species": [ { "name": "B", "initial_count": 1, "site_type": "SIT", "energy_ref": -3.25 } ], "events": [ { "name": "B_SIT_to_SIT", "from_site": "SIT", "to_site": "SIT", "barrier": 2.18, "prefactor": 1e13, "neighbor_distance": 1.0 }, { "name": "B_SIT_to_INT", "from_site": "SIT", "to_site": "INT", "barrier": 3.42, "prefactor": 1e13, "neighbor_distance": 1.0 } ], "simulation": { "total_time": 10.0, "temperature": 900.0, "random_seed": 12345, "output_interval": 1.0, "max_events": 1000000 } }逐字段解释:
"lattice_file":必须与上一步生成的.lat文件名一致;"species":定义掺杂原子。"initial_count": 1表示只放 1 个 B 原子;"site_type": "SIT"强制其起始在替位位点;"energy_ref"是该物种在参考态下的形成能(单位 eV),用于后续热力学平衡判断(如是否析出);"events":定义所有允许的跃迁事件。每个事件需明确from_site和to_site类型(必须在.lat文件中存在)、跃迁势垒barrier(eV)、指前因子prefactor(s⁻¹),以及neighbor_distance(Å)——这是 MMonCa 判断两个位点是否“相邻”的阈值,直接影响邻接矩阵构建;"simulation":控制模拟节奏。"total_time"是目标模拟时长(秒),注意:kMC 时间是真实物理时间,不是步数;"temperature"决定玻尔兹曼因子;"output_interval"指定每间隔多少秒输出一次快照(.xyz格式)。
注意:
barrier和prefactor必须来自 DFT 计算或可靠文献。MMonCa不做任何势垒预测,它只做速率计算与随机采样。若你只有实验拟合的扩散系数 D,可用阿伦尼乌斯公式反推barrier:D = prefactor * a² * exp(-Ea/kT) / 6(a 为晶格常数),其中prefactor通常取 1e13–1e14 s⁻¹ 量级。
2.4 运行与验证:用mmonca执行并检查输出
一切就绪后,执行主程序:
./mmonca --input=input.json --output_dir=./out成功运行后,./out/下会生成:
trajectory.xyz:按output_interval输出的原子位置序列,可用 OVITO 或 VMD 可视化;statistics.log:每步事件类型、发生时间、当前系统能量等摘要;final_state.json:最终构型的 JSON 描述(含每个原子位置与类型)。
验证是否真正启动:打开statistics.log,首行应为# time(s) event_name site_id_from site_id_to energy_change(eV),随后是类似0.000123 B_SIT_to_SIT 42 87 -0.15的记录。若日志为空或报错No valid events found,大概率是neighbor_distance设得过小(导致无邻接边)或barrier为负值(非法输入)。
3. 掺杂扩散模拟的核心配置:三类事件、四组参数与一个不可妥协的约束
MMonCa 的物理保真度,几乎全部压在events数组的定义质量上。它不内置任何经验势函数,所有跃迁行为均由用户显式声明。我们拆解其事件建模的底层逻辑,并给出工业级模拟中必须覆盖的三类事件及其参数设定原则。
3.1 三类必建事件:体扩散、间隙-替位转换、表面吸附/脱附
在半导体工艺仿真中,仅定义“替位→替位”跃迁远远不够。真实掺杂动力学至少包含以下三类事件,缺一不可:
| 事件类型 | 物理意义 | 典型from_site/to_site | 是否必须 |
|---|---|---|---|
| 体扩散跃迁 | 掺杂原子在晶格内部替位位点间的迁移(如 B 在 Si 中沿<100>跳) | "SIT"→"SIT" | ✅ 必须,定义主扩散通道 |
| 间隙-替位转换 | 掺杂原子从间隙态被捕获到替位,或从替位被踢出到间隙(Frank-Turnbull 机制) | "INT"↔"SIT" | ✅ 必须,决定电活性浓度 |
| 表面吸附/脱附 | 外延生长中,气相原子在衬底表面的吸附、表面迁移、脱附(影响形貌) | "SURF"↔"SIT"或"INT" | ⚠️ 若模拟外延,必须启用 |
提示:
"SURF"是特殊 site_type,需在lattice_gen时通过--surface=true参数生成表面位点。默认.lat文件不含 SURF,强行定义 SURF 事件会导致运行时崩溃。
3.2 四组关键参数的物理含义与取值边界
每个事件的四个参数并非独立,它们共同构成阿伦尼乌斯速率Γ = prefactor * exp(-barrier / (k_B * T))。但 MMonCa 对它们有隐含约束:
| 参数 | 单位 | 合理范围 | 超出后果 | 经验建议 |
|---|---|---|---|---|
barrier | eV | ≥ 0.0 | 负值导致 Γ > 1e300,数值溢出崩溃 | 用 DFT 计算,误差控制在 ±0.1 eV 内;若无数据,宁可设为 0.0(表示无障碍)也不瞎猜 |
prefactor | s⁻¹ | 1e12 – 1e14 | <1e10 导致事件几乎永不发生;>1e15 可能掩盖其他慢过程 | 统一取1e13作为起点,后续根据实测扩散系数反调 |
neighbor_distance | Å | ≥ 最近邻距离 × 0.9 | 过小:邻接矩阵稀疏,扩散被阻断;过大:引入虚假长程跃迁 | 对 Si,设1.2(最近邻 2.35 Å,取 0.5 倍);对 GaAs,设1.4 |
energy_ref(物种级) | eV | 依据 DFT 计算的形成能 | 影响系统总能计算,但不影响动力学路径选择 | B 在 Si 中替位:-3.25 eV;P 在 Si 中:-4.12 eV(数据来源:Materials Project) |
特别强调neighbor_distance:它是 MMonCa 构建跃迁图谱的唯一几何判据。程序会遍历所有位点对,若欧氏距离 ≤ 此值,则认为可跃迁,并将该边加入图谱。它不是最大跃迁距离,而是“是否考虑该对位点之间存在跃迁通道”的开关。设为0.5,则只有紧邻位点对被连接;设为5.0,则可能连通次近邻甚至更远——这会极大增加图谱边数,拖慢采样速度,且无物理依据。
3.3 一个不可妥协的约束:所有事件必须满足“位点类型守恒”
MMonCa 在初始化时会对events做静态校验:若某事件的from_site类型在.lat文件中不存在(例如写了"SURF"但.lat无表面位点),程序立即退出并报错Site type 'SURF' not found in lattice。同样,to_site也必须存在。这不是运行时错误,而是编译期校验——意味着你必须先生成含所需位点类型的.lat文件,再写input.json。常见翻车点:想模拟表面吸附,却忘了在lattice_gen中加--surface=true,然后死磕input.json里的"SURF"事件,徒劳无功。
4. 外延生长模拟专项:从台阶流到岛状生长的建模要点与参数陷阱
当 MMonCa 的应用从“体扩散”延伸到“外延生长”,问题复杂度陡增:不再只是原子在固定格点间跳,而是要处理开放边界、气相通量、表面重构、台阶边缘动力学等。MMonCa 通过扩展lattice_gen和events语义来支持,但必须理解其建模假设,否则结果毫无物理意义。
4.1 构建含台阶的表面晶格:lattice_gen的--step与--terrace_width参数
外延生长的核心是台阶流(step-flow)机制。MMonCa 要求你显式构造一个带单原子台阶的表面晶格。使用lattice_gen的台阶参数:
./lattice_gen \ --crystal=diamond \ --a=5.431 \ --output=si_step_6x6x1.lat \ --size=6,6,1 \ --surface=true \ --step=true \ --terrace_width=3 \ --step_height=1参数详解:
--surface=true:启用表面位点生成(SURF类型);--step=true:在表面引入单原子台阶;--terrace_width=3:台阶平台宽度为 3 个原子列(沿某个晶向);--step_height=1:台阶高度为 1 个原子层(即单原子高)。
生成的.lat文件中,SURF位点会被标记为不同子类型:SURF_TERRACE(平台区)、SURF_STEP_EDGE(台阶边缘)、SURF_CORNER(台阶拐角)。MMonCa 允许你在events中针对这些子类型定义不同事件,例如:
{ "name": "Si_gas_to_terrace", "from_site": "GAS", "to_site": "SURF_TERRACE", "barrier": 0.8, "prefactor": 1e10 }, { "name": "Si_terrace_to_edge", "from_site": "SURF_TERRACE", "to_site": "SURF_STEP_EDGE", "barrier": 0.3, "prefactor": 1e12 }注意:
"GAS"是特殊 site_type,代表气相原子。它不占用.lat文件中的格点,但必须在species中声明("name": "Si", "site_type": "GAS"),且initial_count设为 0(气相浓度由prefactor和total_time隐式控制)。
4.2 气相通量建模:用prefactor编码束流强度,而非添加“气相原子数”
新手最大误区:试图在species中设置"initial_count": 1000来模拟高束流。这是错误的——initial_count只对固相位点有效。气相通量必须通过GAS相关事件的prefactor实现。其物理逻辑是:prefactor正比于气相原子撞击表面的频率。例如:
- 分子束外延(MBE)中 Si 束流为 0.1 ML/s(单层每秒),对应
prefactor ≈ 1e10 s⁻¹(对 6×6 表面); - 化学气相沉积(CVD)中 SiH₄ 分压更高,
prefactor可升至1e11–1e12 s⁻¹。
没有通用换算表,必须根据你的反应室几何、束流口径、分压值自行估算。一个实用技巧:先用极低prefactor(如1e8)跑 1 秒,看吸附原子数;再线性放大prefactor,使 1 秒内吸附原子数 ≈ 目标 ML 数 × 表面原子数。
4.3 台阶钉扎与岛状生长:用“吸附能差”驱动形貌演化
真实外延中,原子在台阶边缘的吸附能通常比平台区高 0.2–0.5 eV,导致原子优先在台阶处捕获,实现台阶流生长。若想模拟岛状生长(2D nucleation),则需反转这一能差:让平台区吸附能更高,促使原子在平台成核。这通过energy_ref实现:
"species": [ { "name": "Si", "initial_count": 0, "site_type": "GAS", "energy_ref": 0.0 } ], "events": [ { "name": "Si_gas_to_terrace", "from_site": "GAS", "to_site": "SURF_TERRACE", "barrier": 0.6, "prefactor": 1e11, "energy_ref_change": -3.8 // 吸附后能量下降 3.8 eV }, { "name": "Si_gas_to_edge", "from_site": "GAS", "to_site": "SURF_STEP_EDGE", "barrier": 0.4, "prefactor": 1e11, "energy_ref_change": -4.2 // 吸附后能量下降 4.2 eV → 更稳定 } ]energy_ref_change是事件专属参数,表示该跃迁发生后,系统总能的变化(eV)。它不改变energy_ref,但影响后续事件的热力学可行性(如脱附概率)。正是这个差值(-4.2 vs -3.8)让台阶边缘成为能量洼地,原子自然流向那里。
5. 避坑指南:MMonCa 用户最常踩的 4 个血泪现场与自救方案
MMonCa 功能强大,但因其领域专用性,新手极易陷入一些隐蔽深坑。这些坑往往不报错,或报错信息模糊,导致调试数日无果。以下是我在多个模拟项目中反复验证的 4 个高频致命问题,按“现象→原因→解决”结构列出,句句来自真实翻车记录。
5.1 现象:模拟运行飞快,statistics.log里全是同一类事件,其他事件从不发生
原因:barrier值单位错误或数量级错位。例如,把 DFT 给的 2.18 meV 误当 2.18 eV 输入,导致exp(-0.00218/kT)≈ 1,而其他事件 barrier=2.18 eV 的exp(-2.18/kT)≈ 1e-10,速率差 10 个数量级。
解决:立刻检查所有barrier值是否为eV单位。用 Python 快速验证:import math; print(math.exp(-2.18/(8.617e-5*900)))→ 应得~1e-12;若输0.00218,结果是~0.998。养成习惯:DFT 输出务必标注单位,粘贴前用编辑器全局搜索meV替换为eV。
5.2 现象:mmonca启动后立即报错Segmentation fault (core dumped),无其他日志
原因:lattice_gen生成的.lat文件与input.json中site_type不匹配,且错误发生在内存分配阶段(如请求 0 个SURF_STEP_EDGE位点但代码未做空检查)。
解决:第一步,用head -n 20 si_step_6x6x1.lat查看文件末尾是否有SURF_STEP_EDGE字样;第二步,确认lattice_gen命令中--step=true和--surface=true同时启用;第三步,删掉build/重新 cmake(旧 build 缓存可能残留错误头文件)。90% 的 segfault 源于此。
5.3 现象:trajectory.xyz显示原子在跳,但final_state.json中所有原子 site_id 都是 0,或全为同一值
原因:neighbor_distance设得过大,导致图谱中每个位点都连向所有其他位点,kMC 采样失去局部性,原子随机瞬移到任意位置,违背晶体物理。
解决:将neighbor_distance改为晶体最近邻距离的 0.8–0.9 倍。查 Si 最近邻:2.35 Å,故设1.9;查 GaAs:2.48 Å,设2.0。运行后检查statistics.log中site_id_from和site_id_to的差值,95% 应在 ±50 内(对 6×6×1 表面,总位点约 200,差值 >100 说明乱跳)。
5.4 现象:模拟耗时远超预期,1 秒物理时间跑了 2 小时,CPU 占用 100%
原因:events数组中存在大量低 barrier、高 prefactor 的事件,导致每步采样需遍历数千条边,且多数被拒绝(因随机数未命中)。典型如定义了 50 种SIT→INT事件,但只有一种物理合理。
解决:精简events。用grep -c '"from_site": "SIT"' input.json统计 SIT 出发事件数,若 >5,必有冗余。删除所有 DFT 未验证、文献未报道的事件。保留原则:每个from_site类型,最多 3 种to_site类型(如 SIT→SIT、SIT→INT、SIT→SURF)。MMonCa 不是穷举工具,是精准建模工具。
6. 进阶验证与可信度加固:用三组对照实验锁定模型偏差源
跑出一个能动的模拟只是起点,真正的工程价值在于:这个“能动”的结果,多大程度上反映了真实物理?我在某高校合作的 SiC 掺杂项目中,曾用三组低成本对照实验,在一周内定位出模型偏差根源——不是 DFT 数据不准,而是lattice_gen的对称性处理忽略了 SiC 的双原子基元。以下是我的标准化验证流程,你可直接复用。
6.1 对照组 A:零温度极限测试(T → 0 K)
目的:验证事件 barrier 排序是否与 DFT 结果一致。
方法:将input.json中"temperature"改为1.0(K),"total_time"设为0.001,运行。此时exp(-barrier/kT)差异极大,最低 barrier 事件将垄断所有跃迁。
预期结果:statistics.log中 99% 以上事件应为 barrier 最小的那个。例如,若B_SIT_to_SITbarrier=2.18 eV,B_SIT_to_INTbarrier=3.42 eV,则前者应占绝对主导。
若不符:说明barrier值输错,或.lat文件中位点类型标签错(如把SIT写成SIT1)。
6.2 对照组 B:无限稀释极限测试(single dopant)
目的:排除多体相互作用干扰,聚焦单原子动力学。
方法:确保input.json中"initial_count": 1,且lattice_file尺寸足够大(如8x8x8),使原子间距 >5 nm,杜绝邻居效应。
预期结果:扩散系数D应严格满足阿伦尼乌斯律ln D = ln(prefactor * a² / 6) - barrier / (k_B * T)。用 3 个温度(700K, 900K, 1100K)各跑一次,拟合直线斜率应 ≈-1/k_B(≈ -12000 K)。
若斜率偏差 >5%,检查prefactor是否统一,或a(晶格常数)是否在lattice_gen中输错。
6.3 对照组 C:对称性自洽测试(mirror lattice)
目的:验证晶格生成与事件定义是否尊重晶体对称性。
方法:生成两个镜像.lat文件:si_mirror_x.lat(x 方向镜像)和si_orig.lat,用完全相同的input.json运行。
预期结果:两者的statistics.log中,各事件发生频次分布应完全相同(允许浮点误差)。例如B_SIT_to_SIT在两者中占比均为 82.3±0.1%。
若显著不同(如 82% vs 45%),说明lattice_gen的镜像操作破坏了位点类型标签,或neighbor_distance在镜像后导致邻接关系畸变——此时应回退到lattice_gen源码,检查--mirror参数实现。
我的习惯是:每次新搭一个材料体系(如从 Si 切换到 GaN),必跑这三组对照。它们不产出论文图,但能省下两周调试时间。有一次,对照组 C 揭示出
lattice_gen对 wurtzite 结构的 c/a 轴处理有 bug,我直接给作者提了 PR,三天后 merged。技术债越早还,后面路越宽。希望帮到你。
本文还有配套的精品资源,点击获取