求助，在molpro用多参考态进行结构优化如何加速收敛

TY5272611

大家好，初用Molpro，在molpro用CASPT2做结构优化和能量计算，基组为cc-pvtz和cc-pvdz外推至CBS，活化空间7e6o，一周后仍未收敛，请问有什么解决办法可以加速收敛？构型和输入输出文件如下

zjxitcc

你文字描述说“基组为cc-pvtz和cc-pvdz外推至CBS”，但输入文件似乎是cc-pVTZ和cc-pVQZ外推至CBS，与你所说不完全相符。

另外，您似乎对CASSCF计算需要合理的活性轨道、CASPT2解析导数的计算量没有经验和概念，计算输入十分粗暴，从输出内容看基本是nonsense，像是学习3天的Molpro新手想做3年老司机的计算。我需要问一些基本问题，以防您一开始就走错了：您想做的是中性、二重态的(7e,6o)计算吗？您这个确定需要使用CASPT2做结构优化吗，是审稿人/导师强制要求，还是你自己想做？谁告诉您这个体系需要使用CASPT2做结构优化？您是否对CASSCF/CASPT2计算有一定了解？（我并非嘲笑您，而是需要对这个问题有一定的背景了解，尤其是体系电荷和自旋多重度）

TY5272611

zjxitcc 发表于 2025-8-16 15:58
你文字描述说“基组为cc-pvtz和cc-pvdz外推至CBS”，但输入文件似乎是cc-pVTZ和cc-pVQZ外推至CBS，与你所说 ...

感谢老师的批评，回答您的问题：1.基组cc-pVTZ和cc-pVQZ外推至CBS，是我打错了不好意思；2.是中性二重态的7e6o计算；3.该构型为OH位点H提取的过渡态，T1诊断值为0.07存在多参考特性，采用CASPT2做优化是审稿人的要求；4.对CASSCF/CASPT2计算原理了解还不够深入

zjxitcc

当前(7e,6o)设置不明。建议在CASPT2(5e,5o)/def2TZVP级别下做结构优化，1天内可以算完。活性空间包含
O-H键成键和反键轨道 (2e,2o)
O上的单电子(1e,1o)
O-O过氧键的成键和反键轨道 (2e,2o)
---
合计(5e,5o)。如果你理解并使用下述计算技巧，可以自己把基组换成cc-pVTZ或def2TZVPP，计算耗时可能增加几个小时。

首先我们对初始结构做一个自动CASSCF单点计算，找出重要的轨道，确定活性空间，顺便看看程序自动确定的活性空间是否合适（如果合适，计算就更简单了）。输入文件C5H13O3.gjf内容如下

1. %nprocshared=64
   
2. %mem=160GB
   
3. #p CASSCF/def2TZVP
   
4. 
   
5. mokit{npair=25}
   
6. 
   
7. 0 2
   
8. C     3.038899     -0.629479      0.117225
   
9. H     3.574008     -0.086487     -0.666247
   
10. H     3.381958     -1.665629      0.102518
    
11. H     3.324613     -0.195093      1.078263
    
12. C     1.529433     -0.534959     -0.087424
    
13. H     1.266518     -0.996109     -1.048284
    
14. H     1.009823     -1.093762      0.697877
    
15. C     1.030691      0.909000     -0.083916
    
16. H     1.478088      1.454695     -0.922427
    
17. H     1.352860      1.408979      0.836807
    
18. C    -0.488100      1.014337     -0.167931
    
19. H    -0.863146      0.460077     -1.048545
    
20. C    -0.974147      2.468745     -0.313558
    
21. H    -0.549740      2.901103     -1.221173
    
22. H    -2.061558      2.500671     -0.370583
    
23. H    -0.645583      3.052304      0.548221
    
24. O    -1.110281      0.501376      0.960046
    
25. H    -1.529194     -0.675520      0.798023
    
26. O    -1.586968     -1.739428      0.548966
    
27. O    -1.537549     -1.790067     -0.830260
    
28. H    -0.680916     -2.206150     -1.000839
    
29. 
    

复制代码

核数和内存可以自行修改。npair=25表示除了8个1s芯轨道 对其余25个双占据轨道都构造对应的反键轨道。若直接提交到当前节点，可运行如下Linux命令

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

复制代码

若使用队列脚本提交到计算节点，那就把上面命令写进脚本，并去掉末尾的后台运行&符号。CASSCF计算需要良好的初始轨道，往往不会使用R(O)HF正则轨道作为初始轨道，而MOKIT的automr通过多步计算工作流可自动构造良好的初始轨道，并自动调用量化软件进行CASSCF轨道优化。该任务12 min完成，打开输出文件C5H13O3.out可以看到

1. CASSCF(3e,3o) using program pyscf

复制代码

说明automr自动确定该体系基态最小活性空间为(3e,3o)。再看这两个重要的波函数文件

GVB轨道	C5H13O3_uhf_uno_asrot2gvb25_s.fch
CASSCF自然轨道	C5H13O3_uhf_gvb25_CASSCF_NO.fch

可以用GaussView/Multiwfn打开上述fch文件，观看自然轨道和自然轨道占据数，无需使用Molpro推荐的冷门可视化程序。从xxx_CASSCF_NO.fch文件可以看出(3e,3o)包含O上的单电子 和 O-O过氧键的成键和反键轨道，这部分我们直接接受，认为这3个就是最重要的轨道。由于研究问题是H转移反应，还应当把O-H键的成键和反键轨道放入活性空间，合在一起构成(5e,5o)。xxx_s.fch文件中的轨道是CASSCF初始轨道，我们从这个文件出发，将O-H键（第15号和53号轨道）调换至目标位置，这可通过如下的Python脚本实现

1. from mokit.lib.gaussian import permute_orb
   
2. from shutil import copyfile
   
3. 
   
4. old_fch = 'C5H13O3_uhf_uno_asrot2gvb25_s.fch'
   
5. new_fch = 'C5H13O3_uhf_uno_asrot2gvb25_p.fch'
   
6. copyfile(old_fch, new_fch)
   
7. permute_orb(new_fch, 15, 32)
   
8. permute_orb(new_fch, 53, 36)
   

复制代码

运行python脚本就是简单地

1. python permute.py

复制代码

随即获得文件C5H13O3_uhf_uno_asrot2gvb25_p.fch，内含调换好的轨道，我们读取这套轨道进行CASSCF(5e,5o)单点计算，输入文件C5H13O3_2.gjf内容如下

1. %nprocshared=64
   
2. %mem=160GB
   
3. #p CASSCF(5,5)/def2TZVP
   
4. 
   
5. mokit{ist=5,readno='C5H13O3_uhf_uno_asrot2gvb25_p.fch'}
   

复制代码

提交automr任务的命令不再重复演示。该任务耗时6 min，正常结束后获得文件C5H13O3_uhf_uno_asrot2gvb25_p_CASSCF_NO.fch，内含收敛的CASSCF(5e,5o)轨道。运行命令

1. fch2com C5H13O3_uhf_uno_asrot2gvb25_p_CASSCF_NO.fch

复制代码

随即产生两个文件

Molpro轨道文件	c5h13o3_uhf_uno_asrot2gvb25_p_.a
Molpro输入文件	C5H13O3_uhf_uno_asrot2gvb25_p_CASSCF_NO.com

其中已写好体系坐标、基组数据和读取轨道的关键词。fch2com是MOKIT的一个小程序，用于传轨道给Molpro。我们打开.com输入文件，把末尾3行删掉，改为下述4行CASPT2结构优化关键词

1. save,mo,2140.2,ORBITALS}
   
2. {CASSCF;closed,31;occ,36;CANONICAL;NoExtra}
   
3. {RS2;CORE,8}
   
4. {optg,Gaussian}
   

复制代码

最后一行的Gaussian表示使用量化软件Gaussian的结构优化收敛限。最后提交Molpro计算任务，按你的习惯来。我使用的提交命令见下方，仅供参考

1. molpro -W ./ -t 1 -n 64 -m 300m C5H13O3_uhf_uno_asrot2gvb25_p_CASSCF_NO.com

复制代码

Molpro输出文件为C5H13O3_uhf_uno_asrot2gvb25_p_CASSCF_NO.out，前几分钟可看到类似如下输出

1. Second-order MCSCF: L-BFGS accelerated
   
2. 
   
3. ITER  MIC  NCI  NEG     ENERGY(VAR)    ENERGY(PROJ) ENERGY CHANGE       GRAD(0)  GRAD(ORB)   GRAD(CI)      STEP      TIME
   
4.    1   80   20    0   -421.53606975    -421.53607096   -0.00000121    0.00019530 0.00000022 0.00000000  0.88E-02    109.04
   
5.    2   68   19    0   -421.53607107    -421.53607118   -0.00000012    0.00022699 0.00000134 0.00000000  0.25E-02    225.48
   
6.    3   84   23    0   -421.53607119    -421.53607119   -0.00000000    0.00001963 0.00000000 0.00000000  0.18E-03    344.95
   
7. 
   
8. CONVERGENCE REACHED!  Final gradient:    0.00000001 ( 0.51E-08)
   

复制代码

表明我们传给Molpro的是收敛的CASSCF(5,5o)轨道，所以第一帧结构的CASSCF立即收敛、进入CASPT2能量和力计算，然后经过47步后结构优化收敛。紧接着下方Current geometry是最后一帧结构的坐标，达到目的

1.   47  -423.07806736  -423.07806736    -0.00000000  0.00000579  0.00001079  0.00000143  0.00024086  0.00036354  0.00004815  38136.74
   
2. 
   
3. END OF GEOMETRY OPTIMIZATION.
   
4. 
   
5. Geometry written to block  1 of record 700
   
6. 
   
7. Current geometry (xyz format, in Angstrom)
   
8. 
   
9.    21
   
10. RS2/USERDEF  ENERGY=-423.07806736
    
11. C           2.5832545848       -0.8170285380        0.2241508155
    

复制代码

若还想获得最后一帧结构的CASSCF自然轨道，用于文章中展示或分析，有两种方式：
（1）在Molpro输入文件末尾再加一行{put,xml;keepspherical}；
（2）不改.com文件，直接复制最后一帧结构的坐标，写一个gjf文件，按照上面展示的技巧，先做一个自动CASSCF计算，再调换轨道做CASSCF(5e,5o)计算，获得xxx_CASSCF_NO.fch文件，同样可以用GaussView/Multiwfn打开可视化，进行各种分析。