高斯和ORCA算(N2)2 casscf单点能计算结果不同

隔壁老张sir

高斯的输入文件：
%nprocshared=12
%mem=24GB
%chk=relax_bentTshape_cas12_ddn.chk
# casscf(12,12) maug-cc-pvtz

Title Card Required

0 1
N                   0.47544795    0.00000000   -0.27450000
N                  -0.47544795    0.00000000    0.27450000
N                   0.00000000    0.00000000  -11.25000000
N                   0.00000000    0.00000000   -8.75000000



ORCA的输入文件：
! maug-cc-pv(t+d)z veryTightSCF
%pal nprocs=16 end
%maxcore 1000
%casscf nel 12
                norb 12
end
*xyz 0 1
N                  0.47544795    0.00000000   -0.27450000
N                 -0.47544795    0.00000000    0.27450000
N                  0.00000000    0.00000000  -11.25000000
N                  0.00000000    0.00000000   -8.75000000
*

各位老师，我在orca手册里没找到maug-cc-pvtz，只有maug-cc-pv(t+d)z，这个应该是和高斯中的maug-cc-pvtz是一样的吧？
高斯得到的能量为-217.91564202hartree，ORCA得到的能量为-217.691297156747hartree，为什么这俩软件得到的能量不同呢？

zjxitcc

隔壁老张sir 发表于 2023-4-24 15:55
邹老师您好，还有个问题想请教一下，我看到wjc404回复此贴的最后说“如果真要看两个N2之间的相互作用，用 ...

因为NEVPT2比CASPT2晚提出、晚发展，纯粹是历史原因，这二者间，未来只会是NEVPT2。现在还在用CASPT2的人，要么是明知有NEVPT2却不想换出舒适区，要么是消息过于闭塞不知今年是2023年。

隔壁老张sir

zjxitcc 发表于 2023-4-17 11:52
您最近发的相关帖子均可在我这一层得到回复和解决，我就不一一搬运回复了：
http://bbs.keinsci.com/threa ...

邹老师您好，还有个问题想请教一下，我看到wjc404回复此贴的最后说“如果真要看两个N2之间的相互作用，用CASSCF是不够的，至少NEVPT2。”我想问一下可以用CASPT2吗？因为我参考的这篇文献(Paukku Y, Yang K R, Varga Z, et al. Global ab initio ground-state potential energy surface of N4[J]. The Journal of chemical physics, 2013, 139(4): 044309.)用的是CASPT2，但是CASPT2不是存在一些缺点吗(您在此贴中http://bbs.keinsci.com/forum.php ... 1&fromuid=44565说到过，其他老师也说过)，为什么文章中还会用到这种方法，而不是用CASPT2对标的NEVPT2方法呢？如果我用NEVPT2方法是不是更好呢？

隔壁老张sir

wzkchem5 发表于 2023-4-18 04:43
接近平衡键长的时候MP2可以合理描述。只要学过Moller-Plesset微扰理论，必然知道为什么

好的，谢谢老师

wzkchem5

隔壁老张sir 发表于 2023-4-17 09:07
老师您好，我看您说“N3和N4的距离拉近了以让MP2能够合理描述”，然后将这里改成了“-9.8400 ...

接近平衡键长的时候MP2可以合理描述。只要学过Moller-Plesset微扰理论，必然知道为什么

隔壁老张sir

wjc404 发表于 2023-4-17 14:26
是ORCA的CASSCF收敛到了错误的波函数上。

可以用MP2自然轨道作为CASSCF的初猜轨道。

老师您好，我看您说“N3和N4的距离拉近了以让MP2能够合理描述”， 然后将这里改成了“-9.84000000”，此时的N3和N4的距离是1.09。那如果一开始的坐标就是您改之后的样子，即N                  0.47544795    0.00000000   -0.27450000
N                 -0.47544795    0.00000000    0.27450000
N                  0.00000000    0.00000000  -9.84000000
N                  0.00000000    0.00000000   -8.75000000
那我是不是就不需要改N3的坐标了呢？一般是让N3和N4的距离大致为多少时，可以让MP2能够合理描述呢（大约是1左右吗）？谢谢老师了

隔壁老张sir

wjc404 发表于 2023-4-17 14:26
是ORCA的CASSCF收敛到了错误的波函数上。

可以用MP2自然轨道作为CASSCF的初猜轨道。

十分感谢大佬。计算这个体系主要就是想把这个化学反应过程弄清楚，得到化学反应过程的相关参数

隔壁老张sir

zjxitcc 发表于 2023-4-17 11:52
您最近发的相关帖子均可在我这一层得到回复和解决，我就不一一搬运回复了：
http://bbs.keinsci.com/threa ...

谢谢z神

wjc404

是ORCA的CASSCF收敛到了错误的波函数上。

可以用MP2自然轨道作为CASSCF的初猜轨道。

首先生成MP2自然轨道
1）将下面的代码内容保存为输入文件N_MP2NAT.INP，运行ORCA（N3和N4的距离拉近了以让MP2能够合理描述）。

1. ! RI-MP2 may-cc-pV(T+d)Z aug-cc-pVTZ/C VeryTightSCF PAL4
   
2. %maxcore 2000
   
3. 
   
4. %scf
   
5.   stabperform true
   
6. end
   
7. 
   
8. %mp2
   
9.   natorbs true
   
10. end
    
11. 
    
12. *xyz 0 1
    
13. N                  0.47544795    0.00000000   -0.27450000
    
14. N                 -0.47544795    0.00000000    0.27450000
    
15. N                  0.00000000    0.00000000  -9.84000000
    
16. N                  0.00000000    0.00000000   -8.75000000
    
17. *

复制代码

2）运行以下命令在本地生成N_MP2NAT.mp2nat.molden.input轨道文件

1. cp N_MP2NAT.mp2nat N_MP2NAT.mp2nat.gbw
   
2. orca_2mkl N_MP2NAT.mp2nat -molden

复制代码

3）用molden或multiwfn程序打开N_MP2NAT.mp2nat.molden.input，看最高的6个占据和最低的6个未占据轨道形状是否符合预期。

然后用MP2自然轨道作为CASSCF的初猜，做CASSCF计算
1）将下面的代码内容保存为输入文件N_CAS12.INP，运行ORCA

1. ! may-cc-pV(T+d)Z VeryTightSCF PAL4 MORead
   
2. %moinp "N_MP2NAT.mp2nat"
   
3. %maxcore 2000
   
4. 
   
5. %casscf
   
6.   nel 12
   
7.   norb 12
   
8. end
   
9. 
   
10. *xyz 0 1
    
11. N                  0.47544795    0.00000000   -0.27450000
    
12. N                 -0.47544795    0.00000000    0.27450000
    
13. N                  0.00000000    0.00000000  -11.25000000
    
14. N                  0.00000000    0.00000000   -8.75000000
    
15. *

复制代码

最后得到的能量为-217.915642 a.u.
2）运行以下命令，转化出可以被可视化程序使用的CASSCF轨道文件N_CAS12.molden.input

1. orca_2mkl N_CAS12 -molden

复制代码

3）用molden或multiwfn打开N_CAS12.molden.input查看活性空间的轨道是否符合预期。

以及好奇这个体系的计算意义在哪里。两个N2分子相聚很远，相互作用很弱，可以分别计算。如果真要看两个N2之间的相互作用，用CASSCF是不够的，至少NEVPT2。

zjxitcc

您最近发的相关帖子均可在我这一层得到回复和解决，我就不一一搬运回复了：
http://bbs.keinsci.com/thread-36575-1-1.html
http://bbs.keinsci.com/thread-36314-1-5.html

回复内容：
（1）CASSCF不是黑箱式方法，在常规量子化学程序里往往不能通过一行关键词达到目的，需要多个繁琐步骤。你算这个体系，首先要确定(12,12)里含的是2根N-N sigma单键，4根N-N pi键，不看、不确定就开始算，属于想撞大运。建议阅读点资料，了解多组态方法计算基本常识
(i) 《ORCA CASSCF tutorial》
(ii) 《用Gaussian做CASSCF计算》
(iii) 《CASSCF初始轨道高效构建（一）：局域轨道》
（2）guess=mix是给UHF/UDFT用的，不是给CASSCF用的，不要瞎写。int=acc2e=12是g16默认的，不用写。guess=always和SCF(novaracc,noincfock)此处也用不到。
（3）可以采用fch2mkl小程序从高斯fch文件生成ORCA inp和mkl文件，inp文件内含maug-cc-pvtz基组数据、mkl文件内含轨道，完美解决基组问题、CASSCF收敛问题。
（4）势能面扫描（无论刚性、柔性扫描），强烈建议从 断键 扫至 成键，对你这个体系，就是先从其中一个N...N距离很长的结构，往距离变短的方向扫。
（5）当前扫描坐标较特殊，建议产生一系列坐标文件，挨个算单点，不建议采用高斯自带的势能面扫描。当然，势能面有跳跃 与 GIC功能是否有bug 没有关系。同时，还应该从N...N距离很长的结构，往距离变短的方向算。

下面我用MOKIT做个示例，我们先用automr做个自动的多参考态计算，gjf文件如下

1. %mem=200GB
   
2. %nprocshared=48
   
3. #p CASSCF(12,12)/maug-cc-pvtz guess(fragment=3)
   
4. 
   
5. mokit{}
   
6. 
   
7. 0 1 0 1 0 4 0 -4
   
8. N(fragment=1)   0.51961525    0.00000000   -0.29999998
   
9. N(fragment=1)  -0.51961525    0.00000000    0.29999998
   
10. N(fragment=2)   0.20131400   -0.00002100  -11.35670200
    
11. N(fragment=3)  -0.24553300    0.00002100   -8.64324900

复制代码

提交任务到当前节点（若提交到集群队列，则把下列命令写到提交脚本中）

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

复制代码

耗时仅1 min多一点，得到CASSCF自然轨道文件N2_N2_final_uhf_gvb10_CASSCF_NO.fch。这里我把其中处于平衡位置附近的N2键长稍微拉长了一点，便于让程序确定出的(12,12)就是我们想要的。而且我还用了片段组合波函数作为UHF初猜，保证UHF收敛到了正确的解上、且检验了波函数稳定性。你可以打开*_CASSCF_NO.fch这个文件看自然轨道，会发现automr自动给出的活性轨道就是这个问题里需要的。接着我们算本帖目标结构，gjf文件内容为

1. %mem=300GB
   
2. %nprocshared=48
   
3. %chk=test.chk
   
4. #p CASSCF(12,12)/maug-cc-pvtz nosymm int=nobasistransform guess=read
   
5. 
   
6. title
   
7. 
   
8. 0 1
   
9. N     0.47544795    0.00000000   -0.27450000
   
10. N    -0.47544795    0.00000000    0.27450000
    
11. N     0.00000000    0.00000000  -11.25000000
    
12. N     0.00000000    0.00000000   -8.75000000
    

复制代码

将上述自然轨道fch文件转化为chk文件，即

1. unfchk N2_N2_final_uhf_gvb10_CASSCF_NO.fch test.chk

复制代码

提交g16任务，即

1. g16 test.gjf &

复制代码

这样就会自动从*_CASSCF_NO.fch这个文件中读取收敛的CASSCF轨道，投影到当前结构上。由于两个结构很相近，高斯CASSCF用16圈即可收敛，能量为 -217.915642 a.u. 下面我们在ORCA中重复出一样的结果。运行

1. fch2mkl test.chk

复制代码

此举会产生inp和mkl文件，将mkl转化为gbw文件，供inp文件读取，即

1. orca_2mkl test_o -gbw

复制代码

打开test_o.inp文件，将前几行修改成CASSCF(12,12)的计算关键词，如

1. %pal nprocs 24 end
   
2. %maxcore 13000
   
3. ! VeryTightSCF
   
4. %casscf
   
5. nel 12
   
6. norb 12
   
7. maxiter 200
   
8. ActOrbs NatOrbs
   
9. CI
   
10.   MaxIter 200
    
11. end
    
12. end
    
13. %scf
    
14. Thresh 1e-12
    
15. Tcut 1e-14
    
16. end
    

复制代码

提交ORCA任务，可以看到CASSCF在两三圈后收敛，能量仍-217.915642 a.u. 两个程序一致。另外，如果你有一系列结构要算，我建议用PySCF做CASSCF计算，比上述两个程序会更快一些。用ORCA也还不错，毕竟已经有了很好的初始轨道。