关于CH2CHOH+NO2体系进行NEVPT2计算时CASSCF轨道被转出去的疑问

400__LUX

各位老师好，我在使用多参考态方法做CH2CHOH+NO2体系一个反应时，使用了MOKIT分别对TS和反应物进行NEVPT2/cc-pVQZ的计算。没指定活性空间时，MOKIT分别对TS和反应物选择了（3,3）和（5,5）作为活性空间，随后我将TS的活性空间扩大至（5,5）来保证活性空间大小一致，两者计算的GVB轨道和CASSCF轨道如下文件所示:
链接: https://pan.baidu.com/s/1EmHg4G2el13Y4U3n7_FI0A?pwd=suvr
通过使用Multiwfn观察GVB 轨道，我认为过渡态的15号轨道和22号轨道应该交换，但我在交换轨道后发现，CASSCF计算后的轨道并没有和我预期的一样将15号轨道纳入活性空间。请问各位老师，我应该如何操作来保证两者的活性空间一致？

zjxitcc

在讨论活性空间问题之前，我想问一下你这是什么基元反应的过渡态呢？建议用语言描述一下，或者用ChemDraw之类的软件给出基元反应步骤。如果是N...H...O氢转移过渡态的话，这咋键长看上去不太合理呢？

400__LUX

zjxitcc 发表于 2025-6-8 20:09
在讨论活性空间问题之前，我想问一下你这是什么基元反应的过渡态呢？建议用语言描述一下，或者用ChemDraw之 ...

老师，我这个反应是NO2从CH2CHOH的OH上抽氢的过渡态，即CH2CHOH+NO2→CH2CHO+HNO2。过渡态是在M06-2X-D3/6-311++G**级别上找到的，IRC如下：

zjxitcc

我用TS结构做了freq计算和波函数稳定性检验，在M06-2X-D3/6-311++G(d,p)级别下过渡态结构确实是这样的。这个反应用UCCSD(T)算可能更为简单一些；相比之下，用CASSCF或更高的方法，则难度要大很多，以下是我个人建议的计算步骤

Step 1. 过渡态结构的中等基组单点计算
写一个输入文件，不妨称为ts.gjf，内容如下

1. %nprocshared=32
   
2. %mem=64GB
   
3. #p CASSCF/cc-pVDZ
   
4. 
   
5. mokit{}
   
6. 
   
7. 0 2
   
8. C                  1.83488913    0.85386706   -0.35031703
   
9. H                  1.21731609    1.10500908   -1.20135909
   
10. H                  2.53582118    1.58876512    0.02058300
    
11. C                  1.73135512   -0.37062703    0.26502802
    
12. H                  2.37759917   -0.62272604    1.10537108
    
13. O                  0.86901006   -1.26803709   -0.08208801
    
14. H                 -0.15084101   -0.80922006   -0.20823301
    
15. N                 -1.22152609    0.00497100   -0.00701200
    
16. O                 -2.35271917   -0.37631803   -0.14079301
    
17. O                 -0.86962506    1.11984708    0.32843802

复制代码

提交任务

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

复制代码

算化学反应的话，我通常会先做一些测试，cc-pVDZ基组让计算10 min结束，这时候有UNO轨道（_uno.fch结尾），局域配对UNO轨道（_asrot.fch结尾），GVB轨道（_gvb17_.fch结尾），CASSCF轨道（_CASSCF_NO.fch结尾）。

Step 2. 反应物结构的中等基组单点计算
写一个输入文件，不妨称为r.gjf，内容如下

1. %mem=64GB
   
2. %nprocshared=32
   
3. #p CASSCF/cc-pVDZ
   
4. 
   
5. mokit{npair=17}
   
6. 
   
7. 0 2
   
8. C                 -1.18652274   -0.24849797    0.00000000
   
9. H                 -1.11329937   -1.31598965    0.00000000
   
10. H                 -2.14760933    0.22183459    0.00000000
    
11. C                 -0.06200659    0.50782919    0.00000000
    
12. H                 -0.13522995    1.57532087    0.00000000
    
13. O                  1.22243624   -0.12074610    0.00000000
    
14. H                  1.11201688   -1.07437480    0.00000000
    
15. N                  9.69454932    0.14839092    7.01076901
    
16. O                 10.69403640   -0.08698433    6.42683425
    
17. O                  9.07884138   -0.39318215    7.86103847

复制代码

没指定npair=17的话，在GVB轨道那一步默认只会算出15对轨道，这与过渡态结构的17对GVB轨道不符，不利于我们后面挑选轨道 统一R与TS的活性空间。因此这里手动指定npair=17，automr会增加2对比较合理的轨道。

Step 3. 观察两个结构的轨道，统一R与TS的活性空间
首先肯定是看_CASSCF_NO.fch文件中的自然轨道，发现自动计算出的结果无法实现统一活性空间的目的，接着就是看GVB轨道、局域配对UNO轨道，再结合反应机理（氢转移）来判断。我个人的建议是，对于反应物（以下用符号r或R简称），活性轨道可以取O-H sigma键，C-C pi键，O的孤对电子（平行于C=C-O平面的那个孤对，这个孤对电子轨道后续会与相邻的碳成C=O pi键）。每根键有其成键轨道和反键轨道，因此3根键对应CAS(6,6)，还必须考虑二重态的那个单占据轨道，加起来就是CAS(7,7)。在实际计算中我发现C-O sigma键总是往活性空间里面挤，试图把其他1个轨道给交换出去。索性我们把C-O sigma键也囊括进活性空间，最终取CAS(9,9)。但是我描述的这些轨道的位置并不全在第20~28号位置，需要先调换一下轨道顺序，python脚本如下
from mokit.lib.gaussian import permute_orb
from mokit.lib.rwwfn import read_nif_from_fch, read_na_and_nb_from_fch, write_eigenvalues_to_fch
import numpy as np

fchname= 'r_uhf_uno_asrot.fch'
permute_orb(fchname,18,22)
permute_orb(fchname,30,26)
permute_orb(fchname, 9,21)
permute_orb(fchname,39,27)
permute_orb(fchname, 8,20)
permute_orb(fchname,40,28)
nif = read_nif_from_fch(fchname)
na, nb = read_na_and_nb_from_fch(fchname)
ev = np.zeros(nif)
ev[:nb] = 2.0
ev[nb:na] = 1.0
write_eigenvalues_to_fch(fchname, nif, 'a', ev, True)
（论坛代码框有时候有bug，加了代码框可能会吞掉代码）
这样就把目标轨道全放置在第20~28号位置了。这个过程不会产生新文件，只会更新r_uhf_uno_asrot.fch文件。上述脚本中涉及到的轨道序号是我自己看了r_uhf_uno_asrot.fch文件后得出的，你也需要看一下轨道，才能理解这个脚本的含义。接着我们用这套轨道做CASSCF(9,9)/cc-pVDZ计算，就是写一个r2.gjf文件，内容如下

1. %mem=64GB
   
2. %nprocshared=32
   
3. #p CASSCF(9,9)/cc-pVDZ
   
4. 
   
5. mokit{ist=5,readno='r_uhf_uno_asrot.fch'}

复制代码

在我电脑上计算发现CASSCF花了9圈就收敛了，并且能量下降平缓

1. macro iter   1 ( 21 JK    4 micro), CASSCF E = -357.020313477913  dE = -2.55893026e-02  S^2 = 0.7500000
   
2. macro iter   2 ( 21 JK    4 micro), CASSCF E = -357.028259735942  dE = -7.94625803e-03  S^2 = 0.7500000
   
3. macro iter   3 ( 21 JK    4 micro), CASSCF E = -357.044733120973  dE = -1.64733850e-02  S^2 = 0.7500000
   
4. macro iter   4 ( 21 JK    4 micro), CASSCF E = -357.059339396232  dE = -1.46062753e-02  S^2 = 0.7500000
   
5. macro iter   5 ( 17 JK    4 micro), CASSCF E = -357.060509583794  dE = -1.17018756e-03  S^2 = 0.7500000
   
6. macro iter   6 ( 16 JK    4 micro), CASSCF E = -357.060700579413  dE = -1.90995619e-04  S^2 = 0.7500000
   
7. macro iter   7 ( 14 JK    4 micro), CASSCF E = -357.060735754023  dE = -3.51746101e-05  S^2 = 0.7500000
   
8. macro iter   8 (  7 JK    2 micro), CASSCF E = -357.060738239960  dE = -2.48593739e-06  S^2 = 0.7500000
   
9. macro iter   9 (  4 JK    1 micro), CASSCF E = -357.060738262315  dE = -2.23549819e-08  S^2 = 0.7500000
   

复制代码

说明上述9个活性轨道大概率都保留了下来，没有在轨道优化过程中被旋转出活性空间。对于过渡态（以下用符号ts或TS简称），按照上面说的9个轨道，我们也是要先调换一下轨道顺序，python脚本如下

1. from mokit.lib.gaussian import permute_orb
   
2. permute_orb('ts_uhf_uno_asrot.fch',11,23)
   
3. permute_orb('ts_uhf_uno_asrot.fch',37,25)
   
4. permute_orb('ts_uhf_uno_asrot.fch',16,22)
   
5. permute_orb('ts_uhf_uno_asrot.fch',32,26)

复制代码

然后写ts2.gjf进行CASSCF(9,9)/cc-pVDZ计算

1. %mem=64GB
   
2. %nprocshared=32
   
3. #p CASSCF(9,9)/cc-pVDZ
   
4. 
   
5. mokit{ist=5,readno='ts_uhf_uno_asrot.fch'}

复制代码

也是很快收敛、能量平稳下降。这样我们就在DZ基组下 花费相对比较少的时间 统一了R与TS的活性空间。

Step 4. cc-pVQZ基组下的多参考计算
我们将上述CASSCF(9,9)/cc-pVDZ收敛的轨道投影至CASSCF(9,9)/cc-pVQZ，这一步也用python脚本完成，以R结构的轨道文件为例，脚本内容为
from mokit.lib.gaussian import proj2target_basis
proj2target_basis(fchname='r_uhf_uno_asrot_CASSCF_NO.fch',target_basis='cc-pVQZ',nmo=28)
（论坛代码框有时候有bug，加了代码框可能会吞掉代码）
前28个轨道是双占据+活性轨道，这个轨道数量不同于传统ROHF的占据轨道数量，因此需要专门指明。这个操作会得到r_uhf_uno_asrot_CASSCF_NO_proj.fch文件，内含cc-pVQZ基组下的CASSCF(9,9)初始轨道。接着我们可以做更高等级的计算，例如FIC-NEVPT2，输入文件r3.gjf内容如下

1. %mem=170GB
   
2. %nprocshared=64
   
3. #p NEVPT2(9,9)/cc-pVQZ
   
4. 
   
5. mokit{ist=5,readno='r_uhf_uno_asrot_CASSCF_NO_proj.fch',FIC,NEVPT2_prog=ORCA}

复制代码

这样automr会先调用PySCF做CASSCF(9,9)/cc-pVQZ做轨道优化，3圈左右收敛，然后传轨道给ORCA，进行FIC-NEVPT2计算。对于TS结构的操作是一样的，先用python脚本做轨道投影，然后写ts3.gjf文件读取轨道进行QZ基组下的CASSCF和FIC-NEVPT2计算。

九曜

zjxitcc 发表于 2025-6-19 23:28
我用TS结构做了freq计算和波函数稳定性检验，在M06-2X-D3/6-311++G(d,p)级别下过渡态结构确实是这样的。这 ...

Step 3 这里使用的是居于配对的UNO做后续的操作，使用GVB是否同样可以完成后面的计算呢。是应该同时检查这两个轨道，哪个更“合理”，符合化学直觉，就用哪个做后续的操作么。

zjxitcc

九曜 发表于 2025-6-20 16:24
Step 3 这里使用的是居于配对的UNO做后续的操作，使用GVB是否同样可以完成后面的计算呢。是应该同时检查 ...

一般是这样：如果GVB轨道没有发生sigma-pi轨道混合（产生banana bonds），我就用GVB轨道做后续计算；如果发生了sigma-pi轨道混合，那有两种做法：使用lin_comb_two_mo反向混合banana bonds，得到sigma-pi分离的轨道，然后使用这组轨道；不使用lin_comb_two_mo函数，直接用更前面的局域配对UNO（这个默认是用PM局域化的，不会发生sigma-pi混合）。这里为了尽量让脚本看起来更简单一些，我就没用lin_comb_two_mo函数，用了局域配对UNO。

400__LUX

zjxitcc 发表于 2025-6-20 16:37
一般是这样：如果GVB轨道没有发生sigma-pi轨道混合（产生banana bonds），我就用GVB轨道做后续计算；如果 ...

感谢老师这么细致的解答，按照上文方法自己计算了一下，的确能够选择一致的活性空间。

CCLI

想问楼主为什么不用计算产物？如果要做，是不是也要在同一个活性空间？

400__LUX

CCLI 发表于 2025-6-21 10:29
想问楼主为什么不用计算产物？如果要做，是不是也要在同一个活性空间？

是的，的确得在同样活性空间内计算产物能量，方便后面计算速率常数。

CCLI

zjxitcc 发表于 2025-6-19 23:28
我用TS结构做了freq计算和波函数稳定性检验，在M06-2X-D3/6-311++G(d,p)级别下过渡态结构确实是这样的。这 ...

我有几个问题想问老师：1.基组选择，在MOKIT计算Step1是否得和orca或者高斯得水平相等，例如我是在def-TZVP水平下计算的，gjf也要写这个吗？
2.MOKIT自动选取活性空间，轨道位置是自动调整的吗？
3.楼主为什么只考虑到反应物和过渡态活性空间统一，而不考虑产物(后面楼主回应我确实要考虑)
4.调整轨道的脚本是否适用于一般的化学反应？
5.对于闭壳层单重态和开壳层单重态，用step1的输入写法，MOKIT能够区别开吗？

zjxitcc

CCLI 发表于 2025-6-23 16:54
我有几个问题想问老师：1.基组选择，在MOKIT计算Step1是否得和orca或者高斯得水平相等，例如我是在def-TZ ...

1. def-TZVP基组在gjf文件里写作TZVP，这是Gaussian的规定，MOKIT使用的gjf文件也沿用这一规定。在Step 1中，MOKIT的automr程序会调用Gaussian进行RHF和UHF计算，仍然使用TZVP。
2. 重要的轨道自动放在HONO和LUNO附近，按轨道占据数排序。例如把H2O分子中两根O-H键拉长，xxx_gvb4_s.fch文件中HONO-1和HONO是两根O-H键的成键轨道，LUNO和LUNO+1是两根O-H键的反键轨道，自动判断需要CASSCF(4,4)计算。但在计算化学反应时，如果反应物/过渡态/产物均采用automr独立进行多组态计算，很可能算出不一样的的活性空间。为保持不同结构下使用同一大小的活性空间，经常需要自己调整轨道顺序，加入一些轨道，就如3L中多次使用的permute_orb()那样，调整后读取轨道再次进行多组态计算。
3. 楼主已回答。
4. permute_orb()函数可用于任意反应任意物种的轨道的调换，但被调换的轨道序号 显然是取决于具体化学反应、具体基组的。
5. MOKIT的automr程序自几年前诞生以来 一直都会自动考虑闭壳层单重态与开壳层单重态，这是自动多组态/多参考计算必备基本技能之一。你用automr做一个自动CASSCF计算，看automr输出文件便能看出做了几步什么计算，不是那种“RHF计算->RHF正则轨道当做CAS初始轨道->CASSCF计算”低端操作。

CCLI

zjxitcc 发表于 2025-6-23 23:40
1. def-TZVP基组在gjf文件里写作TZVP，这是Gaussian的规定，MOKIT使用的gjf文件也沿用这一规定。在Step 1 ...

好的，谢谢老师

CCLI

zjxitcc 发表于 2025-6-23 23:40
1. def-TZVP基组在gjf文件里写作TZVP，这是Gaussian的规定，MOKIT使用的gjf文件也沿用这一规定。在Step 1 ...

老师好，我照着您的步骤，根据我的体系，将CASSCF(2,2)/cc-pVDZ收敛的轨道投影至CASSCF(2,2)/def2TZVP，出现如下的错误，我的脚本写法如下
from mokit.lib.gaussian import proj2target_basis
proj2target_basis(fchname='2TIM-3-10-OCH3-ts-bs-CAS-cc-pVDZ_uhf_uno_asrot_CASSCF_NO.fch',target_basis='def2TZVP',nmo=144)


已经解决，麻烦老师了