如何计算四自由基特征值y1（LUNO+1值）？输出的.log文件如何查找“LUNO+1”值?

Cchemistry

新手求助，目前为止了解到两种方法计算：第一种方法，用MOKIT软件直接输出四自由基特征值y1（LUNO+1值）。第二种方法，直接用gaussian计算，在生成的log文件中获得四自由基特征值y1（LUNO+1值）。

下面具体介绍一下所了解的两种四自由基特征值y1（LUNO+1值）计算方法。

第一种方法，参考http://bbs.keinsci.com/thread-49325-1-1.html这篇帖子中zjxitcc老师的回答“双自由基特征y0值就是LUNO占据数，四自由基特征y1值就是LUNO+1占据数。建议使用MOKIT自动做这种计算，默认就会输出双自由基特征y0，四自由基特征y1”。我尝试使用这种方法计算。

使用Gaussview导出.gjf文件，下面是BBTX1a-aaa-c.gjf内容

%nprocshared=28

%mem=28GB

#p CASSCF/cc-pVDZ

mokit{}

0 5

C                 -2.50501200    2.41956700   -0.71437600

C                 -0.78217200    0.89792000   -0.61449900

C                 -3.64131700    3.26564500   -0.80545000

C                 -1.33203200    3.12269000   -0.35540900

C                 -0.25889500    2.19430300   -0.29027800

C                  0.19381200   -0.10147300   -0.58265300

C                 -3.33047000    4.57251500   -0.52210900

H                 -4.64010600    2.93446700   -1.06347800

S                 -1.63716300    4.80915900   -0.13427300

S                  1.42710900    2.16888800    0.05186100

C                  1.47841600    0.38097300   -0.23861000

H                 -3.99245900    5.42724100   -0.51385300

C                  2.67366700   -0.35595500   -0.11826100

C                  2.73748400   -1.79484800   -0.35134400

C                  3.93684000    0.25859700    0.24725600

C                  3.99871100   -2.53716500   -0.21692900

N                  1.71994600   -2.56731400   -0.69476200

C                  5.17851400   -0.50924700    0.37304100

N                  4.09272800    1.54633500    0.49915600

C                  5.21209300   -1.90352300    0.14233600

N                  3.88988000   -3.84319800   -0.46277800

S                  2.29915900   -4.13429100   -0.84234300

N                  6.23709500    0.22511800    0.71714700

S                  5.71071000    1.79473000    0.87577000

H                  6.13123400   -2.46773600    0.23614600

C                 -2.18683900    1.01191900   -0.88427400

C                 -3.10472000   -0.05299600   -1.29593600

H                 -3.20251200   -0.23722800   -2.36572700

C                 -3.87047900   -0.85453700   -0.41542500

C                 -3.81333700   -0.69693800    0.99828900

C                 -4.73492000   -1.85861000   -0.93739000

C                 -4.57299000   -1.49365300    1.84597800

H                 -3.15769800    0.06146800    1.41610500

C                 -5.50119600   -2.66325600   -0.10482300

H                 -4.78835000   -1.99195700   -2.01518300

C                 -5.41361100   -2.47196600    1.28590600

H                 -6.01076900   -3.09831100    1.94549500

C                 -4.50220600   -1.32443400    3.34589500

H                 -5.49062600   -1.11293500    3.76912900

H                 -4.13374400   -2.23587400    3.83033900

H                 -3.83493200   -0.50499100    3.62492000

C                 -6.41303100   -3.72819400   -0.66852800

H                 -6.12668300   -4.72440500   -0.31268600

H                 -7.45238700   -3.56438900   -0.36219600

H                 -6.38401700   -3.74264100   -1.76116600

H                  0.03059700   -1.14908600   -0.78989500

将该文件放入xftp中，并在Xshell中输入命令，具体输入命令如下“ automr BBTX1a-aaa-c.gjf >BBTX1a-aaa-c.out 2>&1”，回车，输出的.log文件显示“Normal”正常，但是搜索关键词“LUNO”，显示“查找不到”，没有找到LUNO+1的值。

第二种方法，考虑到要算四自由基特征值y1，只需要得到自然占据轨道LUNO+1的值即可，就参考了帖子http://bbs.keinsci.com/thread-3529-1-1.html中sobereva老师回答的“如果用高斯，在后HF计算时写pop=no就能给出自然轨道及其占据数。”来算自然轨道占据数。

%nprocshared=28

%mem=28GB

%chk=BBTX1a-aaa-b.chk

#P UB3LYP/6-31G(d,p) pop=no

Title Card Required

0 5

C                 -2.50501200    2.41956700   -0.71437600

C                 -0.78217200    0.89792000   -0.61449900

C                 -3.64131700    3.26564500   -0.80545000

……（同上，忽略）

输出的.log文件显示“Normal”正常，但是搜索关键词“LUNO”，显示“查找不到”，没有找到LUNO+1的值。

求助各位老师，我的操作是否有误，该如何修改？如果想要找到.log文件中的“LUNO+1”的值，应该搜索什么关键词？

hebrewsnabla

automr 输出文件里面有。zou老师说了“默认就会输出四自由基特征y1”。实际上输出的是 tetraradical character(2c^2) y1，无论你搜索 四自由基(tetraradical)，还是特征(character)，还是 y1，都能找到。搜不到的话从头至尾目视输出文件也能看到啊，又不长。

注意你指定的automr输出文件是BBTX1a-aaa-c.out，不是什么.log

wal

1.你要看的是MOKIT的输出，而不是Gaussian计算任务的输出。MOKIT会依次调用很多程序来完成一次多参考计算，我记得默认工作流Gaussian只进行HF计算产生UNO初猜，离拿到多参考波函数还差得远。如果你的MOKIT已经计算完成，则末次计算产生的fch文件中应该保存了CASSCF级别的自然轨道，可以看占据数得到LUNO值，并不是搜索LUNO关键词
2.完全有误，你进行的pop=no拿到的是B3LYP/6-31G(d,p)级别单参考波函数的自然轨道(好像不save甚至都拿不到自然轨道，只做布居分析来着)，与CASSCF多参考波函数毫无关系。。

zjxitcc

补充一点，开壳层单重态物种才有 双自由基特征、四自由基特征。五重态没有这些指标，要么是您计算错了（例如写错了自旋多重度），要么是您误解了双自由基特征、四自由基特征。

Cchemistry

zjxitcc 发表于 2025-6-9 20:39
补充一点，开壳层单重态物种才有 双自由基特征、四自由基特征。五重态没有这些指标，要么是您计算错了（例 ...

非常感谢老师的回复，我没搞清楚四自由基特征的相关知识，不知道算四自由基特征值需要几重态，我把五重态改成单重态后能进行四自由基特征值的计算了
%nprocshared=28
%mem=28GB
#p CASSCF/cc-pVDZ

mokit{}

0 1
C                 -2.50501200    2.41956700   -0.71437600
……（同上，省略）

Cchemistry

hebrewsnabla 发表于 2025-6-9 20:32
automr 输出文件里面有。zou老师说了“默认就会输出四自由基特征y1”。实际上输出的是 tetraradical charac ...

非常感谢老师的回复，根据zjxitcc老师的提醒，能计算四自由基特征值后，我在BBTX1a-aaa-c.out文件中找到了四自由基特征值。
----------------------- Radical index -----------------------
biradical character   (1-2t/(1+t^2)) y0=  0.999
tetraradical character(1-2t/(1+t^2)) y1=  0.542
Yamaguchi's unpaired electrons  (sum_n n(2-n)      ):  7.320
Head-Gordon's unpaired electrons(sum_n min(n,(2-n))):  5.337
Head-Gordon's unpaired electrons(sum_n (n(2-n))^2  ):  4.488
-------------------------------------------------------------

zjxitcc

Cchemistry 发表于 2025-6-12 16:45
非常感谢老师的回复，根据zjxitcc老师的提醒，能计算四自由基特征值后，我在BBTX1a-aaa-c.out文件中找到 ...

注意，每个理论方法都有其对应的y0和y1值，UHF方法有，GVB方法有，CASSCF方法也有。你需要根据你的需求，计算并读取相应的数据。

automr做自动CASSCF计算是分多个步骤的，一个完整的CASSCF计算会依次输出UHF/GVB/CASSCF三个方法下的y0和y1值，你注意不要读错了。什么是读错，我举个例子：比如你想要CASSCF下的y0和y1值，觉得比较严谨，但你没完成automr自动CAS计算，算了一会儿看到UHF的y0值就以为是目标数据了。

Cchemistry

zjxitcc 发表于 2025-6-12 16:53
注意，每个理论方法都有其对应的y0和y1值，UHF方法有，GVB方法有，CASSCF方法也有。你需要根据你的需求， ...

老师，我在MOKIT说明书“https://jeanwsr.gitlab.io/mokit- ... 5.html#4524-fch2mkl”里看到了有关UDFT的计算   
“4.5.24 fch2mkl
Transfer MOs from Gaussian to ORCA. One example is shown below

fch2mkl h2o.fch
This is used for transferring R(O)HF/UHF/CASSCF orbitals. Two files will be generated: h2o_o.inp and h2o_o.mkl. The _o characters are added to avoid file overwritten in case that there is already a xxx.inp file in the current directory (for example, generated by the utility fch2inp). See more examples below:

(1) Common DFT calculations

fch2mkl h2o.fch -dft wB97M-V
fch2mkl h2o.fch -dft 'b3lyp d3bj'
fch2mkl h2o.fch -dft 'HSE06 D3zero'”
老师，我尝试用MOKIT来做UDFT优化并进行四自由基特征值的计算时显示报错，是不是我输入的内容不正确，请老师指正
下面是我的gjf输入文件
“%nprocshared=28
%mem=28GB
#p opt ub3lyp/6-31g(d,p)

mokit{}

0 1
C                  2.73465400    2.92570400    0.00568400
……（同上，忽略）”
下面是报错提示
“ERROR in subroutine parse_keyword: specified method 'ub3lyp' not supported.

All supported methods are GVB, CASCI, CASSCF, DMRGCI, DMRGSCF, NEVPT2,
NEVPT3, CASPT2, CASPT2K, CASPT3, MRMP2, OVBMP2, MRCISD, MRCISDT, MCPDFT,
FICMRCCSD, MkMRCCSD, MkMRCCSD(T), BWMRCCSD, BWMRCCSD(T), BCCC2b, BCCC3b.”

老师，如果不能用MOKIT来做UDFT结构优化和四自由基特征值计算的话，能否像下面这样先用Gaussian进行UDFT（opt ub3lyp/6-31g(d,p)）结构优化，再用MOKIT进行“CASSCF/cc-pVDZ”四自由基特征值计算
下面是Gaussian进行opt ub3lyp/6-31g(d,p)结构优化的输入文件，进行结构优化
“%nprocshared=28
%mem=28GB
%chk=BBTX1a.chk
# opt ub3lyp/6-31g(d,p)

Title Card Required

0 1
C                 -3.54718600   -1.05862871   -3.93953278
……（同上，省略）”
下面是MOKIT进行“CASSCF/cc-pVDZ”四自由基特征值计算的输入文件
“%nprocshared=28
%mem=28GB
#p CASSCF/cc-pVDZ

mokit{}

0 1
C                 -2.50501200    2.41956700   -0.71437600
……（同上，省略）”

hebrewsnabla

Cchemistry 发表于 2025-6-16 20:11
老师，我在MOKIT说明书“https://jeanwsr.gitlab.io/mokit-doc-mdbook/chap4-5.html#4524-fch2mkl”里看 ...

UDFT级别y1: 先用高斯做正确的UDFT计算（注意要考虑对称破缺），然后参考手册的这一节 https://jeanwsr.gitlab.io/mokit- ... c_unpaired_from_fch

EDIT：还是看10L的详细回复吧。

CASSCF级别y1: 饭一口一口吃，先学会UDFT级别再说。

zjxitcc

Cchemistry 发表于 2025-6-16 20:11
老师，我在MOKIT说明书“https://jeanwsr.gitlab.io/mokit-doc-mdbook/chap4-5.html#4524-fch2mkl”里看 ...

很遗憾，您描述的基本上没有对的内容，DFT计算关键词也是错的，缺乏对 开壳层单重态 和 双自由基/四自由基特征 的正确理解，以及对MOKIT功能的认识。在没掌握 开壳层单重态物种的UKS计算 前，就考虑CASSCF计算 是十分危险的事情，学习跳步过大。这里我举一个三角烯的例子，请仔细、反复阅读

Step 1. 单点计算+波函数稳定性检验
Gaussian输入文件triangulene_ini_stab.gjf内容示例如下

1. %chk=triangulene_ini_stab.chk
   
2. %mem=96GB
   
3. %nprocshared=48
   
4. #p UB3LYP/6-311G(d,p) em=GD3BJ nosymm guess=mix scf(xqc,maxcycle=500) stable=opt
   
5. 
   
6. Title
   
7. 
   
8. 0 1
   
9. 三角烯分子的元素和直角坐标
   
10. [空行]
    
11. [空行]

复制代码

提交任务（注意这里是用g16，不是automr）

1. g16 triangulene_ini_stab.gjf

复制代码

这是对三角烯分子的初始结构进行单点计算，获取稳定的UKS波函数。计算级别UB3LYP/6-311G(d,p)只是一个示例，你喜欢用其他泛函、基组也可以，只要是合理的就行。色散校正em=GD3BJ对平面型单分子的结构影响很小，是否考虑该校正 都可以。Gaussian计算完成后，稳定的UKS波函数保存在triangulene_ini_stab.chk文件中，我们复制一份，即运行Linux命令

1. cp triangulene_ini_stab.chk triangulene_ini.chk

复制代码

下一步要用到这个chk文件。

Step 2. 结构优化+频率分析
Gaussian输入文件triangulene_ini.gjf内容示例如下

1. %chk=triangulene_ini.chk
   
2. %mem=96GB
   
3. %nprocshared=48
   
4. #p opt=calcfc freq UB3LYP chkbasis em=GD3BJ nosymm scf(xqc,maxcycle=500) guess=read geom=allcheck
   
5. [空行]
   
6. [空行，这个文件就这么短，不要再写多余的]

复制代码

提交任务

1. g16 triangulene_ini.gjf

复制代码

这样就会从triangulene_ini.chk文件中读取稳定UKS波函数开始结构优化。假设该计算任务顺利完成，接着用GaussView打开triangulene_ini.log文件，确认没有虚频，以此证实算出的是极小点结构（而非鞍点）。接着运行Linux命令复制一份chk文件

1. cp triangulene_ini.chk triangulene_stab.chk

复制代码

下一步要用到riangulene_stab.chk这个文件。

Step 3. 波函数稳定性检验
我们对初始结构检验了波函数稳定性，但结构优化过程中 分子结构改变了，因此我们无法100%保证最终结构——极小点结构 的UKS波函数也是稳定的。保险起见，应当检验一下。Gaussian输入文件triangulene_stab.gjf内容示例如下

1. %chk=triangulene_stab.chk
   
2. %mem=96GB
   
3. %nprocshared=48
   
4. #p UB3LYP chkbasis em=GD3BJ nosymm scf(xqc,maxcycle=500) guess=read geom=allcheck stable=opt
   
5. [空行]
   
6. [空行，这个文件就这么短，不要再写多余的]

复制代码

提交任务

1. g16 triangulene_stab.gjf

复制代码

这样就会读取triangulene_stab.chk文件中的结构和轨道，检验波函数稳定性。计算完成后，运行Linux命令

1. grep 'SCF Done' triangulene_stab.log

复制代码

如果只能搜索出1处SCF Done，并且显示1~3圈内SCF收敛，说明UKS波函数稳定；如果搜索出多处SCF Done，且能量一次比一次低，说明当前UKS波函数不稳定，频率分析（无虚频）结论作废，需要继续优化几何结构。

Step 4. 获取UKS下的y0, y1值
注意，每个理论方法下都有它对应的y0和y1值，你想要哪个就求哪个。这里我们做完了UKS计算，可以顺便求一下它的y0和y1值（如果你不想要，可跳过此步，直接看下一步）。创建一个test.py文件，内容如下

1. from mokit.lib.gaussian import uno
   
2. from mokit.lib.wfn_analysis import calc_unpaired_from_fch
   
3. uhf_fch = 'triangulene_stab.fch'
   
4. uno_fch = 'triangulene_stab_UNO.fch'
   
5. uno(fchname=uhf_fch)
   
6. unpaired_e = calc_unpaired_from_fch(fchname=uno_fch,wfn_type=1,gen_dm=False)

复制代码

提交任务

1. python test.py

复制代码

屏幕上会显示类似如下内容

1. ----------------------- Radical index -----------------------
   
2. biradical character   (1-2t/(1+t^2)) y0=  1.000
   
3. tetraradical character(1-2t/(1+t^2)) y1=  0.022
   
4. Yamaguchi's unpaired electrons  (sum_n n(2-n)      ):  5.302
   
5. Head-Gordon's unpaired electrons(sum_n min(n,(2-n))):  3.760
   
6. Head-Gordon's unpaired electrons(sum_n (n(2-n))^2  ):  2.745
   
7. -------------------------------------------------------------

复制代码

Step 5. CASSCF单点计算，获取相应的y0, y1值
创建一个triangulene.gjf文件，内容如下

1. %mem=170GB
   
2. %nprocshared=64
   
3. #p CASSCF/cc-pVDZ
   
4. 
   
5. mokit{Npair=11}
   
6. 
   
7. 0 1
   
8. 优化好的三角烯分子的直角坐标

复制代码

提交任务（注意，这里不是用g16，而是用automr）

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

复制代码

这样会自动进行CASSCF计算，算完之后在输出文件中看y0和y1。

你不一定真的需要用三角烯作练习，算个萘分子也可以。后续还可以做更多的多组态方法波函数分析，可以阅读MOKIT线上文档Visualization of upaired electrons。如果将来发表文章时使用了MOKIT，记得恰当引用。