关于用三种泛函对自旋转变周期晶体[Fe(pz)][Ni(CN)4]高自旋态做空间结构优化的问题

hahajisuanxuexi

各位老师好：

我希望用CP2K对自旋转变（spin crossover）周期晶体[Fe(pyrazine)][Ni(CN)4]（属于MOF类型材料，晶体结构见图1）高自旋态做空间结构优化。

该材料在常压低温下是低自旋态，为抗磁态，Fe(ii)的d层电子构型是t2g6，自旋为0，Ni自旋为0；常压300K下会自动转变成高自旋态，此时是顺磁态，Fe(ii)的d层电子构型变为t2g4eg2，自旋为2，Ni​自旋为0；空间结构上高低自旋态之间的差别主要体现在Fe-N八面体的键长的不同（高---2.1+A，低---1.9+A），其他化学键键长基本不变，原胞总共20个原子。

自旋转变研究领域有不少研究（比如最近的一个DOI：10.1063/5.0157187）显示tpssh、b3lyp和m06l是三个描述自旋转变材料性质效果较好的泛函，所以我想用这三个泛函分别对[Fe(pyrazine)][Ni(CN)4]的高低自旋态做空间结构优化对比下不同泛函的模拟结果。


最近有研究（10.1021/acs.jpcc.2c01030）对[Fe(pyrazine)][Ni(CN)4]同类型的材料[Fe(pyrazine)][Pt(CN)4]用quantum espresso做了PBE+D2的计算模拟，因此我参照这个研究的方法用QE的PBE+D2对[Fe(pyrazine)][Ni(CN)4]的高低自旋态做了空间结构优化，在对根据实验数据设置的高低自旋态晶体空间结构所对应的输入文件分别设置不同的自旋极化参数后（高自旋输入文件见附件1，低自旋输入文件见附件2），vc-relax计算都成功收敛，并且结果也体现出了高低自旋态的空间结构和磁性质特点，计算出的键长数据符合预期（优化完成的晶胞结构高自旋见图2，低自旋见图3）。

随后我尝试在QE中用tpssh、b3lyp和m06l泛函做空间结构优化。不过对于tpssh和m06l泛函，我没有找到能够使优化正常进行不报错（比如程序运行不到几分钟就会报错“Error in routine cdiaghg (161): S matrix not positive definite”）的赝势组；而b3lyp则出现不论高还是低自旋态虽然SCF收敛但是一直循环“Using ACE for calculation of exact exchange”计算导致超算100+核运行很久都不结束的问题（见附件3）。因此我决定尝试用CP2K进行tpssh、b3lyp和m06l三种泛函的模拟。


我使用Multiwfn生成了高低自旋态空间结构优化的输入文件，其中低自旋态用tpssh（扩胞成2*2*2只算Γ点，附件4）、b3lyp（扩胞成2*2*2只算Γ点）和m06l（用原胞，K点为6*6*6，附件5）三种泛函都运行成功，其中t和b这两个杂化泛函优化结果符合预期（分别见图4、5），m算出来的Fe-N键长虽然处于低自旋态范围，但是本应该相等的键长长度有区别（见图6）。

高自旋态模拟中b3lyp（附件6）虽然能够收敛并且程序运行成功退出，但得到的优化后的结构并不理想，超胞中的8个原胞有一部分处于高自旋态（Fe-N键长大于2），一部分处于低自旋态（Fe-N键长小于2），且同一个原胞中Fe-N键长不全部相等（见图7）。用tpssh优化时一开始我用的设置和b3lyp几乎一致，但是scf不收敛（输出文件中在10E-5左右徘徊），所以将EPS_SCF从10E-6调高到10E-5（附件7）再次尝试计算，这次SCF可以收敛，但是几何优化中的受力（输出文件中的Maximum gradient数据）运行超过12小时都无法收敛到默认的4.5E-4（最接近的时候能到8E-4），查看输出的几何结构变化的数据pdb文件，发现整体结构只在最开始优化的那个阶段有比较明显的变化，到了后期几乎不再改变，而且后期结构出现了使用b3lyp优化时类似的问题（见图8）。用m06l（输入附件8）优化时SCF一般都徘徊在10-2到10-3之间（输出附件9），我尝试了一些基矢和赝势的组合（https://cp2k-basis.pierrebeaujean.net/），比如专用于PBE和MGGA的基矢赝势组，但SCF收敛问题没有解决。

我也尝试了用CP2K的PBE+D3的方法（用的原胞和6*6*6K点）优化高自旋态（附件10），运行10+小时后出现了类似tpssh高自旋态时受力无法达到收敛限的情况，不过输出的几何结构变化轨迹pdb文件到后期几乎处于稳定状态且键长等数据符合高自旋状态（见图9）；另外用2*2*2扩胞的PBE的高自旋态优化结果作为tpssh优化的初猜也会产生到后期受力无法达到收敛限且优化后期结构存在上文所述类似问题（见图10）的情况。

随后我尝试使用CDFT的方式固定Fe的磁矩，但是尝试的两种泛函pbe和tpssh（见附件11）在程序运行一定时间后都会出现“Target value of constraint”和“Current value of constraint”相差悬殊的问题（见附件12）。

请问各位老师，请问我在上文中使用的各种输入文件参数设置是否存在问题？有没有什么方法能够使tpssh、m06l或b3lyp的常规DFT或者CDFT结构优化结果成功收敛到可以令人满意的高自旋态？