计算化学公社

标题: 求助，在molpro用多参考态进行结构优化如何加速收敛 [打印本页]
作者
Author: TY5272611    时间: 2025-8-16 15:09
标题: 求助，在molpro用多参考态进行结构优化如何加速收敛
大家好，初用Molpro，在molpro用CASPT2做结构优化和能量计算，基组为cc-pvtz和cc-pvdz外推至CBS，活化空间7e6o，一周后仍未收敛，请问有什么解决办法可以加速收敛？构型和输入输出文件如下
作者
Author: zjxitcc    时间: 2025-8-16 15:58
本帖最后由 zjxitcc 于 2025-8-16 15:59 编辑

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

另外，您似乎对CASSCF计算需要合理的活性轨道、CASPT2解析导数的计算量没有经验和概念，计算输入十分粗暴，从输出内容看基本是nonsense，像是学习3天的Molpro新手想做3年老司机的计算。我需要问一些基本问题，以防您一开始就走错了：您想做的是中性、二重态的(7e,6o)计算吗？您这个确定需要使用CASPT2做结构优化吗，是审稿人/导师强制要求，还是你自己想做？谁告诉您这个体系需要使用CASPT2做结构优化？您是否对CASSCF/CASPT2计算有一定了解？（我并非嘲笑您，而是需要对这个问题有一定的背景了解，尤其是体系电荷和自旋多重度）
作者
Author: TY5272611    时间: 2025-8-16 22:19
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计算原理了解还不够深入
作者
Author: zjxitcc    时间: 2025-8-19 17:23
本帖最后由 zjxitcc 于 2025-8-19 17:27 编辑

当前(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

  5. mokit{npair=25}

  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

复制代码
核数和内存可以自行修改。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

  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

  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

  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

  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

  3. END OF GEOMETRY OPTIMIZATION.

  5. Geometry written to block  1 of record 700

  7. Current geometry (xyz format, in Angstrom)

  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打开可视化，进行各种分析。





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