优化长程校正泛函w参数的简便工具optDFTw
注：在北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）中专门有一节“DFT泛函的调控”，其中笔者用了50多页幻灯片对DFT泛函的调控做极尽全面系统和深入的讲解，信息量是本文的5倍以上，十分推荐对这部分内容感兴趣的读者参加！
优化长程校正泛函w参数的简便工具optDFTw
optDFTw: A code for optimizing w parameters of long-range corrected functionals
文/Sobereva @北京科音
First release: 2016-Sep-20 Last update: 2026-Oct-30
1 简介
电子间远程作用距离的Hartree-Fock成分为100%而近程距离HF成分为0或较低的泛函被普遍称为长程校正泛函，其中影响结果的关键参数是w（其实是ω，为打字方便就写成w）。针对特定体系优化w参数以获得尽可能好的结果的做法已经被广为使用了，有兴趣者可以看看钟成等人的综述《最优化“调控”区间分离密度泛函理论的研究进展》（DOI: 10.3866/PKU.WHXB201605301）。对w进行优化的一种很好、很普适、被认可度很高的方法是使当前w下计算的HOMO轨道能量尽可能接近电离能的负值，这也被叫做IP-tuning或gap-tuning。这么做的物理思想是对于精确的交换-相关泛函，HOMO能量精确等于电离能的负值，即Koopmans定理严格满足。这么调节w之后，长程校正泛函可以很好地计算激发能、（超）极化率、fundamental gap、单-三重态激发态能量差等问题（但也并非万能，比如有大小一致性问题、JPCB,119,1202发现有的体系的超极化率即便w经过优化还是算不好）。w对体系依赖性大，针对一个体系优化的w，对于另一个体系就往往很不适合，所以对每个被计算的体系原则上都需要优化w，这比普通泛函计算要多花不少代价。
这里介绍笔者开发的极其简单易用的优化w参数的工具optDFTw，以及附带的对w做扫描的工具scanDFTw。它们可以在官网免费下载：http://sobereva.com/soft/optDFTw/。压缩包中exe后缀的是Windows版可执行文件，没后缀的是Linux版可执行文件。压缩包里还附带了示例模板文件和配置文件。
若你的研究中使用这两个程序，请务必引用：Tian Lu, optDFTw program, http://sobereva.com/soft/optDFTw/
这两个程序都支持Gaussian（至少兼容G09和G16）和ORCA（至少兼容6.1.1版）。对任意体系都能用，无论是中性还是带电体系，无论是闭壳层还是开壳层体系。下面依次介绍optDFTw和scanDFTw的使用。
2 optDFTw程序
optDFTw涉及两个量，J和J^2，都是衡量N电子态以及N+1电子态时HOMO能量与电离能的绝对偏差之和。注意J^2不是J简单的取平方
令J或J^2最小化就找到了最优的w。由于J^2对w更敏感，所以optDFTw最小化的是J^2。optDFTw是基于Brent算法来最小化J^2的。优化w的过程是迭代过程，令w收敛到0.0001 Bohr^-1就已经足够精确了，这一般只需要十来圈甚至更少。
使用optDFTw程序前首先要创建一个Gaussian或ORCA的长程校正泛函的单点任务文件作为模板，通常使用基态极小点结构。下面是个典型的Gaussian模板文件，必须叫template.gjf并放在当前目录下
%mem=60GB
%nproc=16
# LC-wPBE/6-311+G**
test
$$
C 0.00000000 0.00000000 -0.52710800
H 0.00000000 0.93885600 -1.11413900
H 0.00000000 -0.93885600 -1.11413900
O 0.00000000 0.00000000 0.67386600
[空行]
[空行]
LC-PBE记得必须写成LC-PBEPBE。LC-PBE0在模板文件里要写成PBEPBE IOp(3/76=1000010000,3/77=1000010000) IOp(3/119=0750000000,3/120=0750000000) IOp(3/130=02500,3/131=02500)。
下面是个典型的ORCA模板文件，必须叫template.inp并放在当前目录下
! LC-PBE def2-SVP tightSCF noautostart miniprint
%maxcore 1000
%pal nprocs 4 end
* xyz $$
O 0.00000000 -0.00001441 -0.34824012
H -0.00000000 0.76001092 -0.93290797
H 0.00000000 -0.75999650 -0.93290797
*
可见模板文件里体系的净电荷和自旋多重度部分必须写成$$，这部分会自动被optDFTw替换。基组可以是混合基组，照常写即可。可以加入帮助SCF收敛的关键词。
optDFTw目录下有个配置文件settings.ini，它应当出现在当前目录下。其中和optDFTw有关的选项如下，应根据需要自行修改
N= 0 1 //N态的电荷和自旋多重度
N+1= -1 2 //N+1态的电荷和自旋多重度
N-1= 1 2 //N-1态的电荷和自旋多重度
idetail= 1 设1会输出每个w值的HOMO轨道能量和体系总能量，设0不输出
icode= 1 设1调用Gaussian，设2调用ORCA
gaupath= g16 调用Gaussian的命令
orcapath= D:\study\ORCA_6.1.1\orca.exe 调用ORCA的命令
wlow= 0.05 w的下限
whigh= 0.7 w的上限
winit= auto w的初值，auto代表取上限和下限的平均值
wtol= 0.0001 w的收敛限
maxit= 100 最大迭代次数
如果optDFTw启动时在当前目录下找不到settings.ini就会使用默认设置，默认设置如上所示。Brent优化算法需要给定w范围以及初猜，给得越合适越可能用越少步数收敛，范围一定要能够把实际的w值囊括在此范围中，默认值通常就是很合适的。
在Windows下启动可以直接双击optDFTw.exe图标，也可以在命令行里启动。Linux下就在命令行中启动，且需要先用chmod命令增加可执行权限。optDFTw启动后会读取当前目录下的settings.ini和模板输入文件，在optDFTw运行过程中会自动产生对应不同电子数的输入文件（Gaussian是N.gjf、N-1.gjf、N+1.gjf，ORCA是N.inp、N-1.inp、N+1.inp），然后从输出文件中读取计算J^2所需的HOMO轨道能量和体系总能量。
optDFTw也允许在命令行方式运行时加上几个选项，格式为：optDFTw [下限] [上限] [初猜] [收敛限]，比如optDFTw 0.2 0.6 0.35 0.0002，这样指定的选项的优先级高于settings.ini中设置的。
下面是实际运行的输出例子，十分易懂，可见每一轮对N、N+1、N-1体系分别算一次单点。经过10轮，最终优化出的w是0.439042 Bohr^-1。
Command of invoking Gaussian: g16
Charge and spin multiplicity of N state: 0 1
Charge and spin multiplicity of N+1 state: -1 2
Charge and spin multiplicity of N-1 state: 1 2
Maximum iterations: 100
Lower limit: 0.050 Upper limit: 0.700 Init w: 0.375 Tol: 0.00010
Initial point:
Running: g16 < N.gjf > N.out
Running: g16 < N-1.gjf > N-1.out
Running: g16 < N+1.gjf > N+1.out
E_HOMO(N): -0.38543 E_HOMO(N+1): 0.03322
E(N): -114.461295 E(N-1): -114.059335 E(N+1): -114.427351
w: 0.375000 J: 0.0172525740 J^2: 0.0002737347
Iteration: 1
Running: g16 < N.gjf > N.out
Running: g16 < N-1.gjf > N-1.out
Running: g16 < N+1.gjf > N+1.out
E_HOMO(N): -0.34906 E_HOMO(N+1): 0.04772
E(N): -114.448114 E(N-1): -114.048881 E(N+1): -114.415805
w: 0.250861 delta-w: -0.124139 J: 0.0655838670 J^2: 0.0027548020
...略
Iteration: 10
Running: g16 < N.gjf > N.out
Running: g16 < N-1.gjf > N-1.out
Running: g16 < N+1.gjf > N+1.out
E_HOMO(N): -0.39990 E_HOMO(N+1): 0.02729
E(N): -114.465217 E(N-1): -114.062189 E(N+1): -114.431092
w: 0.438953 delta-w: -0.000089 J: 0.0099622720 J^2: 0.0000564966
w has converged!
Final w: 0.439042 Bohr^-1 J^2: 0.0000564966
在Gaussian中用长程校正泛函计算时，将IOp 3/107和3/108都设为MMMMM00000就相当于用了w参数为MMMMM/10000的范围分离泛函。Gaussian中使用此例得到的w=0.439042就应当写IOp(3/107=0439000000,3/108=0439000000)了。在ORCA中则应当写%method rangesepmu=0.439042 end。
若optDFTw调用Gaussian或ORCA计算时失败，通常是Gaussian或ORCA计算出错，应自行打开屏幕上提示的它们的out文件，根据末尾的报错信息结合Gaussian或ORCA使用常识自行判断怎么解决。出错>90%的概率是SCF不收敛，参考《解决SCF不收敛问题的方法》（http://sobereva.com/61）试图在模板文件里加帮助收敛的关键词，或者尝试其它基组（优化出的w对基组不是很敏感，用小一点的基组往往更容易收敛）。也可能是净电荷/自旋多重度设置不对、运行环境有问题、内存不够等原因，应自行仔细看报错提示判断。没Gaussian常识不知道怎么解决的话非常建议通过北京科音初级量子化学培训班（http://www.keinsci.com/KEQC）系统学一遍，不熟悉ORCA的话非常建议通过北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）系统学一遍。
3 scanDFTw程序
scanDFTw是在指定范围内按照指定步长对w进行扫描，对每个w会输出J和J^2值。运行前需要以和optDFTw同样的方式编写Gaussian或ORCA模板文件放到当前目录下，并在settings.ini里通过scanlow、scanhigh、scanstep分别设置扫描下限（初值）、上限、步长。也可以用命令行方式运行：scanDFTw [下限] [上限] [步长] [idetail]。比如scanDFTw 0.3 0.5 0.02 1
下面是scanDFTw扫描最后输出的汇总，可见最优的w在0.35附近，即此时J^2是最小的。
w J J2
0.100000 0.1431884760 0.0107253659
0.120000 0.1256572440 0.0083241232
0.140000 0.1094011790 0.0063647484
0.160000 0.0943445150 0.0047806338
0.180000 0.0804025680 0.0035124349
0.200000 0.0675113880 0.0025111313
0.220000 0.0555651750 0.0017318098
0.240000 0.0444868620 0.0011371579
0.260000 0.0341968980 0.0006967073
0.280000 0.0246355040 0.0003846121
0.300000 0.0157222930 0.0001789383
0.320000 0.0082976230 0.0000618425
0.340000 0.0060753980 0.0000185255
0.360000 0.0076470300 0.0000367312
0.380000 0.0144805230 0.0001062586
0.400000 0.0209011120 0.0002185572
0.420000 0.0269540980 0.0003667936
0.440000 0.0326633620 0.0005449820
0.460000 0.0380533580 0.0007479592
0.480000 0.0431441790 0.0009711981
0.500000 0.0479760620 0.0012120297
4 溶剂效应的考虑
为避免一些人犯错，这里特别一提，如果你要在之后的研究中考虑隐式溶剂模型，那么w调控的时候不要直接带着溶剂模型。参考北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）专门深入讲w调控的一节中下面这页ppt。在培训里有具体做法的专门详细的介绍
