怎样计算双自由基特征值？怎样判断是否需要调换轨道？

Cchemistry

求助，读了http://sobereva.com/264这篇博文，按照sob老师所说的方法进行双自由基特征值的计算。我粗浅的理解大致是这样的，但是不知道自己理解的是否正确，请各位老师批评指正。据了解sob老师介绍了两大种方法，

第一种方法是“2 CI波函数下的双自由基特征”，又分两小点，

  第一小点“2.2 CASSCF计算双自由基C4H8实例”用第一步：#P HF/6-31G*第二步：#P CAS(2,2)/6-31G* guess=read方法得到cas的log文件，然后用y=2*c2^2算双自由基特征值。

  第二小点“2.3 用UHF的自然轨道做CAS”用第一步：#P UHF/6-31G* guess=mix第二步：#P CAS(2,2,UNO)/6-31G* guess=read方法得到cas的log文件，然后用y=2*c2^2算双自由基特征值。

第二种方法是“3 UHF波函数的双自由基特征”分三小点，

  第一小点“3.1 基于alpha和beta轨道重叠来计算”，其中

    “(1)基于alpha和beta分子轨道算重叠积分”是用# UHF/6-31G* guess=mix得到的fchk文件获取alpha分子轨道和序号对应的beta分子轨道的重叠积分，取最高的占据的那个轨道的重叠积分作为T计算。

    “(2)基于双正交化轨道算重叠积分”是用# UHF/6-31G* guess=mix得到的fchk文件获得“Singular values of orbital overlap matrix:”的最后一个最小值作为T计算。

  第二小点“3.2 基于<S^2>来计算”是用“#P UHF/6-31G* guess=mix”得到的log文件中最后一个<S^2>进行计算。

  第三小点“3.3 基于UNO占据数来计算”是用#P UHF/6-31G* guess=mix pop=NO得到的log文件中最高占据轨道的Eigenvalues值为T,用y=1-(2*T)/(1+T^2)计算。

其中“2.3 用UHF的自然轨道做CAS”里提到了“这里都是假定与HOMO、LUMO编号对应的那两个UNO是要被纳入活性空间的，但如果和实际情况不符，要取的是别的UNO，那么应该用guess(read,alter)来调换UNO轨道顺序，就像平时调换MO顺序做CASSCF那样。（无论是两步方法还是三步方法，都是这么用alter调换UNO顺序）”，这里面的alter调换轨道是什么意思，是不是用在进行cas计算这一步中的，如何判断所得的双自由基特征值是算的哪一条轨道的，若不是所需轨道该如何调换轨道，如何指定特定的轨道计算其双自由基特征值呢？想请教一下各位老师，劳烦各位老师帮忙看一下，谢谢

sobereva

双自由基特征是对于整体而言的，不是对轨道而言的

guess=alter用来调整轨道顺序，以将需要考虑的轨道纳入活性空间，这是CASSCF最基本计算常识。

如果不会CASSCF的话，北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）有极为全面完整系统详细的讲CASSCF和多参考计算的部分供学习

zjxitcc

推荐使用开源、免费的自动多参考计算程序MOKIT进行CASSCF单点计算，输入文件就是Gaussian的gjf，例如这个C4H8

1. %mem=8GB
   
2. %nprocshared=4
   
3. #p CASSCF/cc-pVDZ
   
4. 
   
5. mokit{}
   
6. 
   
7. 0 1
   
8. C     -0.74742092    1.77656753    0.00000000
   
9. H     -0.62438907    2.32965189    0.92333358
   
10. H     -0.62438907    2.32965189   -0.92333358
    
11. C     -0.74742092    0.30933780    0.00000000
    
12. H     -1.24978876   -0.09122321    0.88294562
    
13. H     -1.24978876   -0.09122321   -0.88294562
    
14. C      0.74742092   -0.30933780    0.00000000
    
15. H      1.24978876    0.09122321   -0.88294562
    
16. H      1.24978876    0.09122321    0.88294562
    
17. C      0.74742092   -1.77656753    0.00000000
    
18. H      0.62438907   -2.32965189   -0.92333358
    
19. H      0.62438907   -2.32965189    0.92333358
    
20. 
    

复制代码

提交任务，就是运行

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

复制代码

体系较小，1 min算完。程序会自动调用Gaussian进行RHF和UHF计算，自动构造活性轨道、自动确定活性空间大小，自动调用PySCF完成CASSCF计算（PySCF是目前世上最强最快CASSCF程序）。这个例子会自动确定为CAS(2,2)。输出文件里有UHF, GVB和CASSCF各种方法下的y0值和未成对电子数，一网打尽，不用自己一个个算、检查

1. E(CASCI)  =      -156.03717089 a.u.
   
2. E(CASSCF) =      -156.03894780 a.u.
   
3. ----------------------- Radical index -----------------------
   
4. biradical character y0 = n_LUNO =  0.729
   
5. tetraradical character y1 = n_{LUNO+1} =  0.000
   
6. Yamaguchi's unpaired electrons  (sum_n n(2-n)      ):  1.853
   
7. Head-Gordon's unpaired electrons(sum_n min(n,(2-n))):  1.457
   
8. Head-Gordon's unpaired electrons(sum_n (n(2-n))^2  ):  1.716
   
9. -------------------------------------------------------------
   

复制代码

完整输出文件可以看附件
test.out (4.68 KB, 下载次数 Times of downloads: 0)

2025-10-17 09:42 上传 Uploaded

点击下载
Click to download

Cchemistry

zjxitcc 发表于 2025-10-17 09:42
推荐使用开源、免费的自动多参考计算程序MOKIT进行CASSCF单点计算，输入文件就是Gaussian的gjf，例如这个C4 ...

老师您好，我用mokit算出来的out文件是这样的，请老师帮忙看看是我哪里算错了。
Enter subroutine get_paired_LMO...
chem_core=4, ecp_core=0, Nskip_UNO=0
Two sets of MOs, ist=1, invoke UNO associated rotation.
$python C4H8-zjxitcc-mokit_uhf_uno_asrot.py >C4H8-zjxitcc-mokit_uhf_uno_asrot.out 2>&1
----------------------- Radical index -----------------------
biradical character   (1-2t/(1+t^2)) y0=  0.757
tetraradical character(1-2t/(1+t^2)) y1=  0.000
Yamaguchi's unpaired electrons  (sum_n n(2-n)      ):  2.027
Head-Gordon's unpaired electrons(sum_n min(n,(2-n))):  1.782
Head-Gordon's unpaired electrons(sum_n (n(2-n))^2  ):  1.940
-------------------------------------------------------------
Leave subroutine get_paired_LMO at Fri Aug 29 21:43:32 2025


Enter subroutine do_gvb...
nbf=96, nif=96, doubly_occ=4, npair=12, nopen=0, nvir=68
LocDocc=F, Lin_dep=F

Warning: AutoMR recommends GVB(12), but user specifies GVB(11). Try to fulfill...
OK, fulfilled. You are supposed to check GVB orbitals after converged.

Updated number of various orbitals:
doubly_occ=5, npair=11, nopen=0, nvir=69

ERROR in subroutine check_gms_path: rungms does not exist.
gms_path=NOT FOUND

zjxitcc

Cchemistry 发表于 2025-10-17 16:28
老师您好，我用mokit算出来的out文件是这样的，请老师帮忙看看是我哪里算错了。
Enter subroutine get_p ...

使用automr做自动CASSCF计算，默认中途会调用GAMESS，需要安装GAMESS，见GAMESS编译教程
MOKIT文档里有汇总常见报错，见https://jeanwsr.gitlab.io/mokit-doc-mdbook/chap_appdx.html