计算化学公社

标题: 高斯和ORCA算(N2)2 casscf单点能计算结果不同 [打印本页]
作者
Author: 隔壁老张sir    时间: 2023-4-17 10:27
标题: 高斯和ORCA算(N2)2 casscf单点能计算结果不同
高斯的输入文件：
%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，为什么这俩软件得到的能量不同呢？

作者
Author: zjxitcc    时间: 2023-4-17 11:52
本帖最后由 zjxitcc 于 2023-4-17 11:54 编辑

您最近发的相关帖子均可在我这一层得到回复和解决，我就不一一搬运回复了：
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)

  5. mokit{}

  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

  6. title

  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也还不错，毕竟已经有了很好的初始轨道。
作者
Author: wjc404    时间: 2023-4-17 14:26
本帖最后由 wjc404 于 2023-4-17 14:33 编辑

是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

  4. %scf
  5.   stabperform true
  6. end

  8. %mp2
  9.   natorbs true
  10. end

  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

  5. %casscf
  6.   nel 12
  7.   norb 12
  8. end

  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。
作者
Author: 隔壁老张sir    时间: 2023-4-17 14:51
zjxitcc 发表于 2023-4-17 11:52
您最近发的相关帖子均可在我这一层得到回复和解决，我就不一一搬运回复了：
http://bbs.keinsci.com/threa ...

谢谢z神
作者
Author: 隔壁老张sir    时间: 2023-4-17 14:52
wjc404 发表于 2023-4-17 14:26
是ORCA的CASSCF收敛到了错误的波函数上。

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

十分感谢大佬。计算这个体系主要就是想把这个化学反应过程弄清楚，得到化学反应过程的相关参数
作者
Author: 隔壁老张sir    时间: 2023-4-17 16:07
wjc404 发表于 2023-4-17 14:26
是ORCA的CASSCF收敛到了错误的波函数上。

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

老师您好，我看您说“N3和N4的距离拉近了以让MP2能够合理描述”， (, 下载次数 Times of downloads: 3) 然后将这里改成了“-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左右吗）？谢谢老师了


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

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

好的，谢谢老师
作者
Author: 隔壁老张sir    时间: 2023-4-24 15:55
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方法是不是更好呢？
作者
Author: zjxitcc    时间: 2023-4-24 16:38
本帖最后由 zjxitcc 于 2023-4-24 16:39 编辑
隔壁老张sir 发表于 2023-4-24 15:55
邹老师您好，还有个问题想请教一下，我看到wjc404回复此贴的最后说“如果真要看两个N2之间的相互作用，用 ...

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



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3