圆锥交叉点(MECI)的优化
圆锥交叉点(MECI)的优化
MECI点是光化学反应中关键的结构，然而S0与S1的CI点由于存在绝热近似失效的问题，难以使用TD-DFT来描述，常常用多参考方法来计算。
CASSCF (Gaussian)
搜索S0/S1交叉点时CAS选项需要写上nroot=2。使用opt=conical时将自动取态平均，因此末尾需要写权重。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
%chk=cas.chk
#p CAS(2,2,nroot=2)/cc-pvdz opt=conical
cas
0 1
C -4.3848464769 -0.1928612943 -0.4063560008
H -4.1069253359 -0.9833660134 0.3125200920
H -5.3273475658 0.0860896888 0.1885275298
C -3.8874570614 1.0629471032 -0.1105968137
H -3.7902077271 1.4961931511 0.9012464564
H -3.6408021428 1.7959965446 -0.8885253636
0.500000000.50000000
收敛时能量应该相当接近：
1
2
3
4
5
6
( 1) EIGENVALUE -77.85594423
( 2) 0.9321043 ( 1)-0.3621449 ( 3) 0.0057068 (
( 2) EIGENVALUE -77.85578500 0.0043 eV
( 1) 0.9319886 ( 2) 0.3621899 ( 3)-0.0146854 (
若遇到大体系，还有一些方法是可以考虑的：
- 纯泛函：对于能量裂分比较大的情况，低HF成分的泛函可以凑合用。但是很难用。
- SF-TDDFT：通过以三重态作参考态避免参套态波函数不稳定的问题，缺点是会引入自旋污染。GAMESS、ORCA、BDF支持。
- MRSF-TDDFT、SA-SF-TDDFT：是前者的升级版，可以在一定程度上解决自旋污染的问题，但MRSF-TDDFT只有GAMESS和OpenQP支持，SA-SF-TDDFT只有Q-Chem支持。
纯泛函(sobMECP)
从原理上讲，DFT没有办法正确描述圆锥交叉点，只能勉强描述圆锥避免交叉。要使用sobMECP搜索圆锥避免交叉，首先要把能量收敛判据给删掉。找到：
1
2
3
4
5
IF (PConv(1) .and. PConv(2) .and. PConv(3) .and. PConv(4) .and. PConv(5)) THEN
Conv = 1
ELSE
Nstep = Nstep + 1
END IF
改为
1
2
3
4
5
IF (PConv(1) .and. PConv(2) .and. PConv(3) .and. PConv(4)) THEN
Conv = 1
ELSE
Nstep = Nstep + 1
END IF
然后准备头文件。对于静态相关比较强的圆锥避免交叉，最好用纯泛函来描述：
1
2
3
4
5
6
7
8
%mem=192GB
%nproc=64
%chk=S1.chk
#n MN15L/6-31g* td force guess(read)
Second State
1 1
由于使用TDDFT计算state2，需要一个extract_energy2文件：
1
2
3
4
5
{
if ($3 == "E(TD-HF/TD-DFT)") {
print $5
}
}
完成后，依次运行：
1
2
3
./prepare.sh
./runfirst.sh
./runMECP.sh
这个不太实用，因为CI点附近势必出现前线轨道简并的情况，而这又将引起Gaussian的L801报错：
1
2
3
**** Fatal Problem: The smallest alpha delta epsilon is 0.57953183D-03
Error termination via Lnk1e in /work/home/gaus12/apprepo/gaussian/16-hy/scripts/../app/g16/l801.exe at Tue Feb 25 10:39:30 2025.
采用IOp(8/11=1)可以避免出现该报错时终止任务，但轨道能量差在计算时需要承担分母，接近0会产生导数数值不稳定的问题，很可能最终无法收敛。
另外，Harvey算法搜索MECI还是时常出问题的，不如BPUPD和惩罚函数等算法稳定。笔者参考XMECP的实现对sobMECP做了一个该版，令其可以实现BPUPD算法的MECI搜索，如有兴趣尝试请前往baneMECP介绍页。
SF-TDDFT
准备好xyz文件，首先进行SF-TDDFT单点计算，确认单重态与三重态的分布：
1
2
3
4
5
6
7
8
9
10
11
! PBE0 def2-SVP RIJCOSX def2/J def2-SVP/C miniprint
%maxcore 3000
%pal
nproc 32
end
%tddft
sf true
nroots 5
end
*xyzfile 0 3 probe.xyz # Spin-Flip计算是从三重态出发的
找到输出部分：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
(SPIN-FLIP GROUND STATE)
STATE 1: E= -0.044023 au -1.198 eV -9661.9 cm**-1 <S**2> = 0.029232 Mult 1
77a -> 76b : 0.992646 (c= 0.99631634)
STATE 2: E= 0.022816 au 0.621 eV 5007.4 cm**-1 <S**2> = 2.011503 Mult 3
76a -> 76b : 0.578826 (c= 0.76080590)
77a -> 77b : 0.412696 (c= 0.64241412)
STATE 3: E= 0.043212 au 1.176 eV 9483.9 cm**-1 <S**2> = 0.096754 Mult 1
76a -> 76b : 0.406651 (c= -0.63769226)
77a -> 77b : 0.563702 (c= 0.75080067)
STATE 4: E= 0.061844 au 1.683 eV 13573.1 cm**-1 <S**2> = 0.996834 Mult 3
74a -> 76b : 0.043170 (c= -0.20777378)
74a -> 77b : 0.020989 (c= 0.14487557)
75a -> 76b : 0.922136 (c= -0.96027893)
STATE 5: E= 0.069465 au 1.890 eV 15245.9 cm**-1 <S**2> = 1.015826 Mult 3
74a -> 76b : 0.916157 (c= -0.95716113)
75a -> 76b : 0.041733 (c= 0.20428661)
75a -> 77b : 0.029991 (c= 0.17317948)
这里STATE 1是SF基态，第一激发单重态是STATE 3，于是使用如下输入文件：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
! CI-opt PBE0 def2-SVP RIJCOSX def2/J def2-SVP/C miniprint
%maxcore 3000
%pal
nproc 32
end
%tddft
sf true
nroots 5
iroot 3 # follow的第一个根
jroot 1 # follow的第二个根
followiroot true # 开启追踪，避免态交叉时follow错
# IROOTMULT TRIPLET # 可以使用此flag来搜索三重态CI
# TDA false # 当前计算单重态CI，可以考虑关掉TDA
end
%geom
maxstep 0.05 # 减小步长
trust -0.05 # 减小步长
end
*xyzfile 0 3 probe.xyz
该方法不是纯态方法，存在自旋污染问题，而且CIopt走到一半后经常污染大到难以分辨哪个是单重态，只能凑合用，比如：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
---------------------
SF-TDA EXCITED STATES
---------------------
the weight of the individual excitations are printed if larger than 1.0e-02
(SPIN-FLIP GROUND STATE)
STATE 1: E= 0.027946 au 0.760 eV 6133.4 cm**-1 <S**2> = 0.147846
96a -> 100b : 0.014731 (c= -0.12137004)
99a -> 100b : 0.044605 (c= 0.21119915)
101a -> 100b : 0.872210 (c= 0.93392172)
101a -> 102b : 0.036161 (c= -0.19016046)
101a -> 104b : 0.014738 (c= 0.12139931)
STATE 2: E= 0.027972 au 0.761 eV 6139.2 cm**-1 <S**2> = 1.169462
99a -> 101b : 0.026155 (c= 0.16172449)
101a -> 101b : 0.921675 (c= 0.96003905)
101a -> 106b : 0.011819 (c= 0.10871489)
可以看到STATE 2的自旋污染很大，<S**2>达到了1.17，指认单重态或是三重态都感觉不太合适，此时搜S0/S1的CI要不要follow这个STATE 2就很难讲了。
MRSF-TDDFT
GAMESS的收敛性不太好，参考态波函数最好由Gaussian准备。首先，要在RODFT下计算一个三重态作为参考态。
1
2
3
4
5
6
7
8
9
10
11
12
%chk=SRR.chk
%mem=96GB
%nprocshared=32
#p ROBHandHLYP/6-31g* nosymm int=nobasistransform
title
1 3
C -0.59562400 -0.17774100 0.00000000
C 0.04669800 -0.09509800 1.24125300
...
这里用的基组一会还要做opt，不要太大
将收敛的波函数传给GAMESS。这一步要用到MOKIT，需要提前装好。
1
fch2inp name.fchk -mrsf
产生的inp文件中，MWords是每个核分配的内存，类似于ORCA的%maxcore，但单位为MW(1MW=8MB)，可能需要修改：
1
2
3
4
5
6
7
8
9
10
11
$CONTRL SCFTYP=ROHF RUNTYP=ENERGY ICHARG=1 MULT=3 NOSYM=1 ICUT=11
MAXIT=200 DFTTYP=BHHLYP TDDFT=MRSF $END
$SYSTEM MWORDS=600 $END
$SCF DIRSCF=.T. $END
$GUESS GUESS=MOREAD NORB=640 $END
$DFT NRAD0=99 NLEB0=590 NRAD=99 NLEB=590 $END
$TDDFT NSTATE=5 $END
$DATA
GAMESS inp file produced by MOKIT,na=136,nb=134,nif=640,nbf=640
...
运行命令：
1
rungms SRR-CI.inp 00 64 > SRR-CI.gms
其中第一个数字是版本号，第二个数字是并行核数。该文件是单点计算，可能耗时稍高于普通TDDFT单点。由于使用了Gaussian传递轨道，SCF通常会在几圈内收敛。然而，由于GAMESS的SCF收敛性并不出色，也可能出现这种情况：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
ITER EX TOTAL ENERGY E CHANGE DENSITY CHANGE DIIS ERROR INTEGRALS SKIPPED
* * * INITIATING DIIS PROCEDURE * * *
1 0 -152.3589223780 -152.3589223780 0.489769120 0.000006459 733451 2859
2 1 -152.0030790890 0.3558432890 0.000001826 0.240699543 721111 4153
3 2 -152.0030821795 -0.0000030905 0.000004500 0.240698185 481673 12927
4 3 -152.0030794241 0.0000027553 0.000003923 0.240699427 502387 12353
5 4 -152.0030777525 0.0000016717 0.000000477 0.240700176 491055 13159
6 5 -152.0030778365 -0.0000000841 0.000001669 0.240700118 430326 16998
7 6 -152.0030784905 -0.0000006540 0.000002828 0.240699829 481864 13880
8 7 -152.0030781570 0.0000003335 0.000000688 0.240700008 496998 13560
9 8 -152.0030780739 0.0000000831 0.000000187 0.240700024 452182 15541
10 9 -152.0030782026 -0.0000001287 1.222392920 0.240700004 414481 16843 < ????
11 10 -151.7930010913 0.2100771113 0.277703111 0.329696390 732215 2890
12 11 -152.0514983457 -0.2584972543 0.442918745 0.228161220 726459 3182
13 12 -152.2967634574 -0.2452651117 0.383817598 0.077095923 729026 3220
14 13 -152.3259835465 -0.0292200892 0.195819508 0.062213954 729549 3251
15 14 -152.3439786445 -0.0179950979 0.297480095 0.050933696 725051 3522
笔者将该GAMESS计算45分钟迭代200圈后收敛失败的结构交给Gaussian，被Gaussian用1分45秒走了30圈拿下，只能说Gaussian不愧是主流量化程序。这种情况下，也许可以尝试切换SCF算法：
1
$SCF DIRSCF=.T. DIIS=.F. SOSCF=.T. $END
注意：GAMESS输入文件的”$”顶头空格必须有，否则不能识别。
如果还是不行，建议求助大师。这里假设顺利运行完毕，在输出文件中找到：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
----------------------------
SUMMARY OF MRSF-DFT RESULTS
----------------------------
STATE ENERGY EXCITATION <S^2> TRANSITION DIPOLE, A.U. OSCILLATOR
HARTREE EV X Y Z STRENGTH
1 NEGATIVE ROOT(S) FOUND.
1 A -1608.7798225490 -1.153 0.0000 0.0000 0.0000 0.0000 0.0000
0 A -1608.7374375718 0.000 (REFERENCE STATE)
2 A -1608.6970167917 1.100 0.0000 0.0000 -0.0000 6.4967 2.3300
3 A -1608.6757530133 1.679 0.0000 -0.1573 -0.0340 0.0000 0.0018
4 A -1608.6528024556 2.303 0.0000 -0.0000 -0.0000 -0.9923 0.0834
5 A -1608.6417893495 2.603 0.0000 -1.0486 -0.1358 0.0000 0.1029
要优化S0与S1的势能面交叉点，即搜索态1和态2的圆锥交叉，更改输入文件，将RUNTYP设为CONICAL，并添加$CONICL字段：
1
2
3
4
5
6
7
8
9
10
11
12
13
$CONTRL SCFTYP=ROHF RUNTYP=CONICAL ICHARG=1 MULT=3 NOSYM=1 ICUT=11
MAXIT=200 DFTTYP=BHHLYP TDDFT=MRSF $END
$SYSTEM MWORDS=600 $END
$SCF DIRSCF=.T. $END
$GUESS GUESS=MOREAD NORB=640 $END
$DFT NRAD0=99 NLEB0=590 NRAD=99 NLEB=590 $END
$CONICL OPTTYP=PENALTY IXROOT(1)=1,2 SIGMA=8.0 $END
$STATPT NSTEP=100 $END
$TDDFT NSTATE=5 $END
$DATA
GAMESS inp file produced by MOKIT,na=136,nb=134,nif=640,nbf=640
...
- NRAD，NLEB，NRAD0，NLEB0均为DFT格点设置，为了降低耗时可以去掉，若SCF不收敛应当加回来再试。
- 如果要加色散校正，可以在$DFT里加IDCVER=4；
- OPTTYP设置的是优化算法，PENALTY是对CI搜索比较(相对)好用的惩罚约束算法。
- SIGMA是惩罚项的拉格朗日乘数，默认为3.5
- IXROOT(1)设置的是要搜索的根。
- NSTEP是步数上限，默认50，一般不够用
- 可以在$STATPT里使用OPTTOL控制几何优化收敛限(默认0.0001)
设置好后，再次提交计算。注意，无论几何优化是否在预期步数内收敛，GAMESS均会进行最终的布居分析并提示ddikick.x: exited gracefully.退出，因此必须手动确认GAMESS是否真的找到了圆锥交叉点：
1
2
3
4
5
6
7
8
9
***** EQUILIBRIUM GEOMETRY LOCATED *****
COORDINATES OF ALL ATOMS ARE (ANGS)
ATOM CHARGE X Y Z
------------------------------------------------------------
C 6.0 -4.3848464769 -0.1928612943 -0.4063560008
H 1.0 -4.1069253359 -0.9833660134 0.3125200920
H 1.0 -5.3273475658 0.0860896888 0.1885275298
...
而不是由于圈数达到上限自动结束几何计算。这玩意感觉很难收敛，也不知道是GAMESS的问题还是圆锥交叉本身难找的问题。可能的解决办法：
- 如果你期望的结构是圆锥避免交叉，去临时文件目录找到输入文件同名的.dat文件，其中有当前结构坐标。把该坐标复制出来，写成RODFT输入文件，重复上述步骤，并调小SIGMA值再提交再次提交计算。
- 如果你期望的结构是圆锥交叉，调大SIGMA值再次提交计算。
- 使用OPTTYP=BPUPD来切换为默认几何算法。注意，此时SIGMA设置无效！
- 从一个较好的初猜结构开始优化。
最后，提供一个MRSF-TDDFT收敛情况绘图的脚本，需要gnuplot环境：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
#!/bin/bash
# 脚本：GAMESS数据处理工具集
# 用法:
# ./ci.sh geom - 提取所有*.gms文件的几何搜索点数据
# ./ci.sh plot - 为所有*.gms文件绘制NSERCH优化过程图
# 显示用法
show_usage() {
echo "用法: $0 <命令>"
echo ""
echo "命令:"
echo " geom 提取所有*.gms文件的几何搜索点数据"
echo " plot 为所有*.gms文件绘制NSERCH优化过程图"
echo ""
echo "示例:"
echo " $0 geom"
echo " $0 plot"
exit 1
}
# 检查参数
if [ $# -lt 1 ]; then
show_usage
fi
COMMAND="$1"
# 几何数据提取功能
geom_function() {
# 检查是否有.gms文件
if ! ls *.gms >/dev/null 2>&1; then
echo "错误: 当前目录下没有找到.gms文件"
exit 1
fi
# 获取所有.gms文件
gms_files=(*.gms)
total_files=${#gms_files[@]}
echo "发现 $total_files 个.gms文件，开始依次处理..."
echo ""
# 依次处理每个.gms文件
for i in "${!gms_files[@]}"; do
LOG="${gms_files[$i]}"
OUTPUT_FILE="${LOG%.gms}_game.xyz"
echo "[$((i+1))/$total_files] 正在处理文件: $LOG"
# 一次性处理，精确复制原脚本逻辑
awk '
BEGIN {
output_file = "'$OUTPUT_FILE'"
# 清空输出文件
system(">" output_file)
atom_count = 0
search_point_count = 0
nserch_count = 0
# 存储所有行内容
delete all_lines
delete search_point_lines
delete nserch_lines
line_num = 0
}
# 存储所有行
{
all_lines[NR] = $0
line_num = NR
}
# 提取原子数目
/TOTAL NUMBER OF ATOMS/ && atom_count == 0 {
match($0, /=\s*([0-9]+)/, arr)
if (arr[1]) {
atom_count = arr[1]
print " 检测到原子数目: " atom_count
}
}
# 记录所有NSERCH行
/NSERCH:/ {
nserch_lines[nserch_count] = $0
nserch_count++
}
# 记录所有几何搜索点行号
/BEGINNING GEOMETRY SEARCH POINT NSERCH/ {
search_point_lines[search_point_count] = NR
search_point_count++
}
# 函数：格式化NSERCH行
function format_nserch_line(line) {
# 解析NSERCH行: NSERCH: 0 E= -876.6553651869 GRAD. MAX= 0.1367237 R.M.S.= 0.0230244
# 提取E值
if (match(line, /E=\s*([-0-9.]+)/, e_arr)) {
e_value = e_arr[1]
} else {
e_value = "N/A"
}
# 提取GRAD. MAX值
if (match(line, /GRAD\.\s*MAX=\s*([-0-9.]+)/, grad_arr)) {
grad_max = grad_arr[1]
} else {
grad_max = "N/A"
}
# 提取R.M.S.值
if (match(line, /R\.M\.S\.=\s*([-0-9.]+)/, rms_arr)) {
rms_value = rms_arr[1]
} else {
rms_value = "N/A"
}
# 返回格式化的字符串
return "MaxF=" grad_max " RMSF=" rms_value " E=" e_value
}
END {
print " 找到 " search_point_count " 个几何搜索点"
# 精确复制原脚本逻辑：只处理偶数编号的点
even_count = 0
for (i = 0; i < search_point_count; i++) {
point_num = i + 1 # 点编号从1开始
# 只处理偶数编号的点
if (point_num % 2 == 0) {
start_line = search_point_lines[i]
# 计算对应的NSERCH行索引
nserch_index = int(i / 2)
if (nserch_index < nserch_count) {
current_nserch = nserch_lines[nserch_index]
# 格式化NSERCH行
formatted_nserch = format_nserch_line(current_nserch)
} else {
formatted_nserch = "MaxF=N/A RMSF=N/A E=N/A"
}
# 写入原子数目和格式化的NSERCH行
print atom_count >> output_file
print formatted_nserch >> output_file
# 根据原脚本，原子数据从第5行开始
data_start = start_line + 5
# 提取原子数据
atoms_found = 0
for (j = 0; j < atom_count && atoms_found < atom_count; j++) {
line_index = data_start + j
if (line_index <= line_num) {
line_content = all_lines[line_index]
# 检查是否为有效的原子数据行
if (match(line_content, /^[[:space:]]*[A-Za-z]{1,2}[[:space:]]+[0-9.-]+[[:space:]]+[-0-9.]+[[:space:]]+[-0-9.]+[[:space:]]+[-0-9.]+/)) {
# 分割字段并格式化输出（移除CHARGE列）
split(line_content, fields)
if (length(fields) >= 5) {
printf "%-2s %15.10f %15.10f %15.10f\n", fields[1], fields[3], fields[4], fields[5] >> output_file
atoms_found++
}
}
}
}
even_count++
}
}
close(output_file)
print " 提取了 " even_count " 个偶数编号的点"
print " 输出文件：" output_file
}
' "$LOG"
echo " ✓ 完成处理: $LOG"
echo ""
done
echo "所有文件处理完成！共处理了 $total_files 个文件。"
}
# 绘图功能
plot_function() {
# 检查是否有.gms文件
if ! ls *.gms >/dev/null 2>&1; then
echo "错误: 当前目录下没有找到.gms文件"
exit 1
fi
# 获取所有.gms文件
gms_files=(*.gms)
total_files=${#gms_files[@]}
echo "发现 $total_files 个.gms文件，开始依次绘图..."
echo ""
# 依次处理每个.gms文件
for i in "${!gms_files[@]}"; do
LOG="${gms_files[$i]}"
FILE_BASE="${LOG%.gms}"
OUTPUT="${FILE_BASE}_nserch.png"
DATA_FILE="${FILE_BASE}_nserch_data.dat"
PLOT_SCRIPT="${FILE_BASE}_plot.gp"
echo "[$((i+1))/$total_files] 正在处理文件: $LOG"
# 提取OPTTOL值
OPTTOL=$(grep "OPTTOL=" "$LOG" | head -1 | sed -n 's/.*OPTTOL=\([0-9.eE+-]*\).*/\1/p')
if [ -z "$OPTTOL" ]; then
echo " 未找到OPTTOL值，使用默认值 0.0001"
OPTTOL="0.0001"
else
echo " 检测到收敛限 OPTTOL=$OPTTOL"
fi
# 首先执行grep命令查看1->2转换
echo " 检查1->2转换:"
grep "1 -> 2" "$LOG" 2>/dev/null || echo " 未找到1->2转换"
echo " 正在提取NSERCH数据..."
# 提取数据
grep "NSERCH:" "$LOG" | awk '{
if(match($0, /NSERCH:[ \t]*([0-9]+)/, step_arr)) {
step = step_arr[1]
}
if(match($0, /E=[ \t]*([+-]?[0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?)/, energy_arr)) {
energy = energy_arr[1]
}
if(match($0, /GRAD\. MAX=[ \t]*([+-]?[0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?)/, grad_arr)) {
grad_max = grad_arr[1]
}
if(match($0, /R\.M\.S\.=[ \t]*([+-]?[0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?)/, rms_arr)) {
rms = rms_arr[1]
}
if(step != "" && energy != "" && grad_max != "" && rms != "") {
print step, energy, grad_max, rms
}
}' | sort -n > "$DATA_FILE"
# 检查数据
if [ ! -s "$DATA_FILE" ]; then
echo " ❌ 未能提取到有效数据，跳过 $LOG"
rm -f "$DATA_FILE"
echo ""
continue
fi
num_points=$(wc -l < "$DATA_FILE")
echo " 数据点数量: $num_points"
# 生成双Y轴gnuplot脚本
cat > "$PLOT_SCRIPT" << EOF
set terminal pngcairo enhanced color font 'Arial,14' size 1200,800
set output '$OUTPUT'
# 设置图表标题和样式
set title "$LOG ciopt (steps: $num_points)" font ',16'
set xlabel "迭代步数 (NSERCH)" font ',14'
set grid
# 设置左Y轴 (对数坐标，梯度和RMS)
set ylabel "梯度值 (对数坐标)" font ',14' textcolor rgb 'red'
set logscale y
set ytics textcolor rgb 'red'
# 设置右Y轴 (线性坐标，能量)
set y2label "能量值" font ',14' textcolor rgb 'blue'
set y2tics
unset logscale y2
# 设置颜色和样式
set style line 1 lc rgb 'red' lw 2 pt 7 ps 1.0 # 梯度最大值
set style line 2 lc rgb 'green' lw 2 pt 9 ps 1.0 # RMS
set style line 3 lc rgb 'blue' lw 3 pt 11 ps 1.2 # 能量值
# 设置图例位置
set key top right
# 绘制图形
plot '$DATA_FILE' using 1:3 with linespoints ls 1 title "GRAD MAX (左轴)" axes x1y1, \\
'$DATA_FILE' using 1:4 with linespoints ls 2 title "R.M.S. (左轴)" axes x1y1, \\
'$DATA_FILE' using 1:2 with linespoints ls 3 title "能量值 (右轴)" axes x1y2, \\
$OPTTOL with lines lw 2 dt 2 lc rgb 'orange' title "收敛限 " axes x1y1
EOF
# 执行gnuplot
gnuplot "$PLOT_SCRIPT"
if [ $? -eq 0 ]; then
echo ""
tail -1 "$DATA_FILE" | awk '{
printf " 最终能量: %.6f\n", $2
printf " 最终梯度: %.2e\n", $3
printf " 最终RMS: %.2e\n", $4
}'
echo " ✓ 绘图完成！输出文件: $OUTPUT"
else
echo " ❌ gnuplot执行失败"
fi
# 清理临时文件
rm -f "$DATA_FILE" "$PLOT_SCRIPT"
echo ""
done
echo "所有文件绘图完成！共处理了 $total_files 个文件。"
}
# 主程序逻辑
case $COMMAND in
"geom")
geom_function
;;
"plot")
plot_function
;;
*)
echo "错误: 未知命令 '$COMMAND'"
show_usage
;;
esac
This post is licensed under CC BY 4.0 by the author.