不同基组和泛函计算得到的含Fe的Fe-C团簇能量排序不一致，如何确定计算级别和结构？

R-S

  最近在计算含过渡金属Fe的Fe-C负离子团簇，选用的是纯泛函PBE和def2-TZVPD，同时在相同基组下用BPW91和TPSS进行了计算对比，结果发现PBE和BPW91得到的几个异构体能量顺序一致，但是TPSS的结果几个异构体的能量顺序发生了变化，即最稳定的异构体由A变为了C。用这三个泛函结合更大的aug-cc-pVTZ基组进行计算发现，几个异构体的能量顺序与def2-TZVPD基组的顺序又不一致，所以现在不知道该用哪个泛函和基组计算的结果来进行分析了。
  请各位帮忙分析一下吧，主要有下面几个问题：
1. 这个FeC负离子团簇是不是用PBE、BPW91、TPSS泛函进行计算是否合理？有哪个泛函是不合适的吗？或者最合适的是哪个？
2. PBE和BPW91比较类似都考虑了电子密度的梯度，TPSS在此基础上又考虑了动能密度，所以TPSS的结果会有明显不同，那么对该体系来说考虑动能密度会更加合理吗？或者说tpss的结果更加合理吗？
3. 相同泛函用aug-cc-pVTZ基组，与def2基组的相对顺序又发生了变化，aug-cc-pVTZ基组更大，所以他的结果更加精确可信吗？ 因为是开壳层，单电子非常多，用CCSD(T)不太好计算，所以用了这种方式，请问大家这样是否存在问题？
4.是否说明这几个异构体能量差别较小，或者说用普通的DFT 很难计算准确，团簇静态相关较强，需要用多参考的方法才能合理描述？
5.所以现在我下一步应该怎么做呢？如何确定用哪个泛函基组和其最稳定异构体结构呢还是该换用其他合理的 方法进行计算呢？

niobium

这里问题很多的，自旋多重度还有波函数是否稳定，不同的泛函会有比较大的差别。这个可以单独做篇文章了。

zjxitcc

你已经是大人了，不用看小学生掰手腕（什么哪个泛函咋咋咋）。如果体系整体带负电荷，基组先定为ma-TZVP；如果还能确定负电荷集中在某个原子或基团上，可以仅给这个原子或基团使用带弥散函数的基组，其他原子不带弥散函数。做个UHF单点计算，检验波函数稳定性，确保波函数稳定，然后看UNO轨道占据数的分布情况，确定体系多组态/多参考特征是否明显。

这里以[FeC]-负离子为例，假设我们想算准它的二重态（其余自旋多重度以此类推），先用不同初猜做几个UHF计算，找到可能最低的、且波函数稳定的UHF解，gjf文件内容较多，以附件上传
FeC_uhf.gjf (3.07 KB, 下载次数 Times of downloads: 10)

2025-4-16 15:56 上传 Uploaded

点击下载
Click to download

FeC_uhf_frag.gjf (3.27 KB, 下载次数 Times of downloads: 8)

2025-4-16 15:56 上传 Uploaded

点击下载
Click to download

提交Gaussian任务，几分钟完成。经过比较，发现E(UHF) = -1300.19156792 a.u., <S^2>= 3.8127可能是能量最低且波函数稳定的UHF解。现在将FeC.chk文件转化为FeC.fch文件，即运行命令

1. formchk FeC_uhf.chk FeC_uhf.fch

复制代码

接着启动python，运行

1. from mokit.lib.gaussian import uno
   
2. uno('FeC_uhf.fch')

复制代码

此举将产生FeC_uhf_UNO.fch文件，用GaussView或Multiwfn打开，第17号轨道为单占据轨道，是1个Fe的3d单电子，符合想象；往占据数减少的方向看，占据数分别是0.785, 0.785, 0.720, 0.116，剩下的都是小于0.02的。这说明在UHF/ma-TZVP下体系表现出明显的六自由基特征，即有极其明显的多参考特征。如果看轨道形状，会发现除了1个单电子，还存在3根Fe-C键。Fe的4个单电子与C的3个单电子呈反铁磁耦合，至少需要CAS(7,7)才能定性合理描述。注意请不要用UKS轨道操作并向我提问。接着我们进行更严谨的CASSCF单点计算，gjf文件如下

1. %mem=150GB
   
2. %nprocshared=64
   
3. #p CASSCF/gen
   
4. 
   
5. mokit{ist=1,readuhf='FeC_uhf.fch'}

复制代码

提交任务，即运行

1. automr FeC.gjf >FeC.out 2>&1

复制代码

在64核、150GB下，这个小体系耗时仅40s，在输出文件中可以看到UHF未成对电子数为6.435，约6个未成对电子，符合上述分析。在40s内发生的事为：automr自动调用GAMESS进行GVB计算，根据GVB自然轨道占据数选出重要的轨道作为CASSCF初始轨道，自动确定活性空间为CAS(9,9)，自动调用PySCF进行CASSCF轨道优化，算完后产生FeC_uhf_gvb10_CASSCF_NO.fch文件，其中含有CASSCF自然轨道及其轨道占据数。用GaussView或Multiwfn打开该fch文件，发现运气好，活性轨道仍然是1个Fe 3d单占据轨道，3根Fe-C键（包括3个成键轨道，3个反键轨道，化学意义清晰），Fe 3d孤对（包括1个3d双占据轨道和1个4d空轨道），加起来等于9个活性轨道。第18~21号自然轨道的占据数为0.600, 0.590, 0.470, 0.095，体现了该体系有中等强度的六自由基特征（与UHF结论稍有不同），难以使用单参考方法（如RHF-based MP2, CCSD(T)）合理描述。若使用能量最低且波函数稳定的UHF/UKS解，双杂化泛函或UCCSD(T)可能可以获得合理结果，但自旋污染会很大，结果需与多参考方法进行比较。如果想囊括更多活性轨道，可以基于GVB轨道采用更大的活性空间进行计算。

以上仅是较为严谨的定性计算、波函数分析。为获得较为准确的电子能量，需用MC-PDFT/CASPT2/NEVPT2/MRCISD等多参考方法计算能量，gjf文件为

1. %mem=150GB
   
2. %nprocshared=64
   
3. #p NEVPT2/gen
   
4. 
   
5. mokit{ist=5,readno='FeC_uhf_gvb10_CASSCF_NO.fch'}

复制代码

将会自动进行CASSCF(9,9)-NEVPT2计算。可以看到，每一次计算都是自动化的，不需要学习Molpro/ORCA/OpenMolcas/PySCF等软件便可轻易调用它们，不需要从零开始，而是读取上一步收敛的轨道进行更高级方法的计算。提交计算任务

1. automr FeC_2.gjf >FeC_2.out 2>&1

复制代码

完成后打开FeC_2.out文件，可以看到E(SC-NEVPT2) = -1300.77130486 a.u. 如果还有另一个结构要算，依葫芦画瓢获得合理的CASSCF(9,9)结果（两个结构需要保证活性空间大小相等）。然后进行NEVPT2计算，获得二者在NEVPT2/ma-TZVP级别下的相对能量。有了这些数据，再去看小学组里谁略胜一筹就比较有说服力了。

R-S

zjxitcc 发表于 2025-4-16 15:56
你已经是大人了，不用看小学生掰手腕（什么哪个泛函咋咋咋）。如果体系整体带负电荷，基组先定为ma-TZVP； ...

非常感谢您这么详细、系统地解答，学习到了很好的处理思路，我准备照您说的思路操作一下，后续有问题的话还希望能继续请教您！

wjc404

原子数量多少？
这种团簇在实际条件下很难孤立存在，容易和外面的原子成键。计算它们的目的是啥？

R-S

wjc404 发表于 2025-4-16 16:22
原子数量多少？
这种团簇在实际条件下很难孤立存在，容易和外面的原子成键。计算它们的目的是啥？

原子数目不多只有十个原子，我们就是主要研究团簇结构和一些电子性质

wjc404

R-S 发表于 2025-4-16 16:26
原子数目不多只有十个原子，我们就是主要研究团簇结构和一些电子性质

普通的NEVPT2和MRCI算不动了，不过DMRG-NEVPT2没准可以。

R-S

wjc404 发表于 2025-4-16 16:41
普通的NEVPT2和MRCI算不动了，不过DMRG-NEVPT2没准可以。

NEVPT2十个原子都很难算的动吗只有四个Fe原子

zjxitcc

R-S 发表于 2025-4-16 16:55
NEVPT2十个原子都很难算的动吗只有四个Fe原子

多组态/多参考方法计算量的主要决定因素是活性空间大小，原子数和基组只是次要因素，这点与单参考方法有很大不同。如果只有1~2个Fe原子，其他都是CHON之类的原子，那NEVPT2或DLPNO-NEVPT2毫无问题。但如果有4个Fe原子，并且它们之间有Fe-Fe键或反铁磁耦合，那么活性空间将超过(20,20)，超过CASSCF方法极限（非堆积机器所能解决），进而导致所有基于CASSCF的MC-PDFT, CASPT2, NEVPT2, MRCISD方法全算不动，只能转而采用DMRG-CASSCF方法。当然了，MOKIT也支持DMRG-CASSCF和DMRG-NEVPT2计算，而且是全自动的（一旦自动确定活性空间大于(16,16)便会自动切换为DMRG计算），但每一步的计算难度和计算量比上面展示的例子大得多，比如UHF那一步可能你就没法算出来一直报错，或者找不到比较低的UHF解。[Fe4S4]之类的大活性空间强关联体系一般是专业多年做多参考计算的课题组才去碰。