钙钛矿分子动力学模拟：力场参数化与相变分析的实战经验

钙钛矿材料（ABO₃型）的相变行为是凝聚态物理里的经典问题。BaTiO₃从立方相到四方相的铁电相变发生在393K，SrTiO₃的反铁畸变相变在105K，PbTiO₃的四方-立方相变在763K——这些相变温度的精确预测是分子动力学模拟的核心挑战。我做钙钛矿分子动力学模拟做了十多年，踩过最多的坑就是力场参数不靠谱导致相变温度偏差200K以上。

有一年我帮一个课题组算BaTiO₃的铁电相变，用的经典Buckingham势参数，跑出来的相变温度是510K，实验值393K，偏差将近30%。查了半天发现是Ti-O势参数的A、B、C三组系数来源不一致——A和B取自一篇1989年的文献，C取自另一篇1996年的文献，两组参数的拟合基准不同，混在一起用必然出问题。重新用同一组拟合参数后相变温度降到405K，偏差控制在3%以内。这个教训让我养成了一个习惯：力场参数必须追溯到一个统一的拟合来源，绝不混用不同文献的参数。

钙钛矿晶体结构与相变序列：建模的物理基础

钙钛矿ABO₃的理想结构是立方钙钛矿（Pm-3m），A位阳离子位于立方体顶点，B位阳离子位于体心，O原子位于面心。实际钙钛矿大多存在结构畸变：BaTiO₃在室温下是四方相（P4mm），Ti原子沿c轴偏移约0.15Å产生铁电极化；SrTiO₃在105K以下发生反铁畸变相变，氧八面体绕[001]轴交替旋转。

钙钛矿分子动力学模拟的第一步是根据研究目标确定要模拟的相变类型。铁电相变涉及B位离子的位移，需要力场能正确描述B-O短程相互作用的各向异性。反铁畸变相变涉及氧八面体的旋转，需要力场能正确描述O-O排斥和A-O相互作用。两种相变的力场要求不同，选错力场类型相变就模拟不出来。

超胞大小的选择取决于相变类型。铁电相变只需要B位离子位移，3×3×3的超胞（135个原子）通常够用。反铁畸变相变涉及氧八面体的交替旋转，旋转周期是晶胞的2倍，至少需要2×2×2的超胞才能容纳一个完整的旋转周期。如果研究两种相变的耦合效应，需要更大的超胞（4×4×4以上）。

力场类型与参数选择：从Buckingham势到壳模型

钙钛矿分子动力学模拟最常用的力场类型是Buckingham势加库仑势。Buckingham势形式是V(r) = A·exp(-r/ρ) – C/r⁶，前一项描述短程排斥，后一项描述色散吸引。配合部分电荷的库仑相互作用，构成Born模型力场。

力场参数的来源有两个途径。一是从文献中引用已拟合好的参数集。BaTiO₃常用的是Tinte等人的参数集（J. Phys.: Condens. Matter, 2003），这个参数集在BaTiO₃的整个相变序列（立方→四方→正交→三方）上都做了验证，相变温度偏差在20K以内。SrTiO₃常用的是Thomas等人的参数集，对反铁畸变相变描述较好。

二是自己做第一性原理拟合。用DFT计算一系列不同构型的能量和力，用GULP或alamode做力场拟合。这个方法的优点是参数和研究体系严格匹配，缺点是工作量大，拟合一个可靠的力场至少需要2-3周。

壳模型（shell model）是Buckingham势的升级版。标准Born模型把原子当作点电荷，极化效应只能通过色散项近似。壳模型把每个原子分成一个带正电的核心和一个带负电的壳，核心和壳之间用弹簧连接，壳的位移模拟原子的电子极化。壳模型对钙钛矿的铁电相变描述明显优于Born模型——BaTiO₃的壳模型力场能正确给出三方→正交→四方→立方的完整相变序列，而Born模型经常只能给出一个相变。

在钙钛矿分子动力学模拟的力场选择中，有一个经验法则：研究铁电性质用壳模型，研究热膨胀和弹性性质用Born模型就够了。壳模型的计算成本是Born模型的2-3倍（因为多了一倍的自由度），对需要长时间模拟的体系，这个成本差异很显著。

LAMMPS中的钙钛矿模拟配置：从data文件到系综选择

LAMMPS是钙钛矿分子动力学模拟的主力工具。建模流程是先用GULP或packmol生成初始结构，导出为LAMMPS的data文件格式，然后在LAMMPS中读入data文件做模拟。

data文件的关键部分是力场参数的写入。Buckingham势在LAMMPS中对应bond_style不对——Buckingham是非键相互作用，用pair_style buck。命令格式是pair_style buck/coul/long 10.0，10.0是截断半径。每对原子类型的参数用pair_coeff指定：pair_coeff 1 2 A rho C，其中1和2是原子类型编号，A、rho、C是Buckingham参数。

壳模型在LAMMPS中的实现稍复杂。LAMMPS没有原生的壳模型pair_style，需要用fix_adapt或用户自定义的壳模型包。另一个选择是用GULP做壳模型计算，GULP原生支持壳模型，且力场参数库更丰富。我的做法是壳模型计算用GULP，常规Born模型计算用LAMMPS。

系综选择取决于研究目标。相变温度测定用NPT系综（恒温恒压），温度从低温逐步升到高温，监控体积和极化的突变点。离子迁移研究用NVT系综（恒温恒体积），在目标温度下做长时间模拟，计算均方位移（MSD）和扩散系数。NPT系综需要选择控压方法——Parrinello-Rahman方法允许晶胞形状变化，适合研究相变引起的对称性变化；Berendsen方法保持晶胞形状不变，适合各向同性热膨胀研究。

时间步长对钙钛矿模拟很关键。含氢原子的体系（如有机-无机杂化钙钛矿CH₃NH₃PbI₃）需要0.5fs的时间步长来正确积分C-H键振动。全无机钙钛矿用1fs足够。如果用壳模型，壳的质量很小，需要更小的时间步长（0.2-0.5fs），否则壳的位置积分会发散。

相变温度的判定方法：序参量与结构因子

钙钛矿分子动力学模拟的核心输出是相变温度的确定。相变温度不是直接读出来的，而是通过监控序参量随温度的变化来判定。

铁电相变的序参量是B位离子的偏移量。在BaTiO₃中，Ti原子沿[001]方向的平均位移〈z_Ti〉在立方相为零，在四方相不为零。做NPT模拟时，每50K或100K一个温度点，每个温度模拟100-200ps达到平衡后再取100-200ps做统计平均。画出〈z_Ti〉随温度的变化曲线，相变温度对应序参量从非零跳变到零的温度。

反铁畸变相变的序参量是氧八面体的旋转角。旋转角通过O原子相对B位原子的偏移计算，在反铁畸变相中相邻八面体的旋转方向相反，旋转角不为零。旋转角的计算用结构因子分析更准确：计算旋转模式对应的波矢q=(1/2,1/2,1/2)的结构因子，在立方相为零，在反铁畸变相不为零。

相变温度判定的精度取决于温度间隔和统计时间。温度间隔太大（>100K）可能跳过相变点，太小（<20K）计算量增加但增益有限。统计时间太短（<50ps）涨落噪声掩盖相变信号。我的标准配置是温度间隔50K，每个温度点平衡100ps+统计100ps，总模拟时间约4-5ns，在64核服务器上跑约48小时。

在钙钛矿分子动力学模拟的项目中，相变模拟有一个容易忽略的问题：超胞尺寸对相变温度的有限尺寸效应。3×3×3超胞给出的相变温度可能比6×6×6超胞偏高30-50K，因为小超胞限制了涨落的相关长度。发表结果前做一个超胞尺寸收敛测试是必要的。

有机-无机杂化钙钛矿的MD模拟：力场选择与稳定性问题

有机-无机杂化钙钛矿（如CH₃NH₃PbI₃，MAPbI₃）的分子动力学模拟比全无机钙钛矿复杂得多，因为有机阳离子（MA⁺或FA⁺）的内部自由度需要正确描述。

MAPbI₃的力场通常采用分层策略：Pb-I框架用Buckingham势或改进的Born-Mayer势，MA⁺分子内部用Harmonic势描述C-N、C-H、N-H键的振动，MA⁺和Pb-I框架之间用Lennard-Jones势加库仑势描述。MA⁺的旋转动力学对材料的光电性质有重要影响——MA⁺在高温四方相中可以自由旋转，在低温正交相中被锁定，这个旋转转变的模拟需要力场能正确描述MA⁺和I⁻框架之间的非键相互作用。

MAPbI₃的MD模拟有一个特殊的稳定性问题：MA⁺中的H原子和I⁻之间的距离在某些构型下可能小于LJ势的排斥墙高度，导致能量发散。解决方法是减小时间步长到0.5fs，或者对H-I距离做约束。另一个问题是Pb-I框架在高温下可能发生非物理的结构崩塌——力场参数对Pb-I排斥势的描述不够硬，需要用排斥更强的势参数或加入角度约束。

钙钛矿MD模拟的经验总结：力场验证是最关键的一步

回过头看，钙钛矿分子动力学模拟的可靠性完全取决于力场质量。不管模拟流程多么规范，力场参数不给力，结果就是错的。我的力场验证标准有三个：第一，室温下的晶格常数和实验值偏差小于2%；第二，已知相变温度和模拟值偏差小于30K；第三，弹性常数和DFT值偏差小于15%。三个条件都满足才用这个力场做生产模拟。

力场验证不通过时有两个选择。一是换一套参数集，不同研究组拟合的参数可能差异很大。二是自己做第一性原理力场拟合，虽然费时间，但对特殊体系（如掺杂钙钛矿、高压相钙钛矿）可能是唯一可靠的选择。

钙钛矿分子动力学模拟的另一个经验是平衡时间要够长。钙钛矿的相变是弱一级相变，相变点附近涨落很大，系统需要很长时间才能达到平衡。我做过一个BaTiO₃的相变模拟，平衡时间从50ps增加到200ps后，相变温度的统计涨落从±30K降到±10K。对于有机-无机杂化钙钛矿，平衡时间需要更长，因为有机阳离子的取向弛豫很慢——MAPbI₃在室温下MA⁺的取向弛豫时间约10-50ps，平衡时间至少200-500ps。

模拟结果的分析也需要注意统计方法的选取。MSD的计算要用Einstein关系式而非Green-Kubo关系式，因为钙钛矿中离子的扩散是跳跃式的，Green-Kubo对跳跃扩散的收敛性不好。极化的计算用Berry相位方法（DFT层面）或偶极矩求和（经典MD层面），两者给出的极化值在壳模型力场下一致性较好。