最近WYD同学做了一个结构优化练习,得到的结构有一个很大的虚频。体系结构如下:
结构优化的输入文件如下:
%nprocshared=36
%mem=64GB
%chk=test.chk
#p opt freq ub3lyp/6-31+g(d) em=gd3bj
Title Card Required
-2 2
B 0.00000000 0.00000000 0.62868400
O 0.21679800 1.28166000 1.47588700
O -1.22781700 0.29637700 -0.27244500
C 0.00000000 2.24101000 0.44449900
O 0.47366000 3.40669300 0.42450900
C -0.89119700 1.63316700 -0.63408300
O -1.26705400 2.21942300 -1.68222700
O 1.22781700 -0.29637700 -0.27244500
O -0.21679800 -1.28166000 1.47588700
C 0.89119700 -1.63316700 -0.63408300
O 1.26705400 -2.21942300 -1.68222700
C 0.00000000 -2.24101000 0.44449900
O -0.47366000 -3.40669300 0.42450900
得到的虚频如下:
其对应的振动模式如下:
振动对应的化学键有两处,振动使得一个键变长,另一个键变短,所以理应可以通过调节两处C-C键的长度来重新构建初始猜测。这样的情况在对称性比较高的结构中常常遇到。本文我们尝试另一种方法:由于对当前任务已经做了频率计算,因此可以从已得到的结构出发重新进行结构优化,并从chk文件读取精确的Hessian矩阵,其输入文件如下:
%nprocshared=36
%mem=64GB
%chk=test.chk
#p opt(rcfc) freq ub3lyp/6-31+g(d) em=gd3bj
geom=allcheck guess=read
这样优化后,顺利得到了无虚频的结构:
得到的结构中两个C-C键的键长不同,也印证了上面的猜想。感兴趣的读者可以尝试通过手动改变键长的方法来重新构建初始猜测。
当然,这种做法也存在着一定的“冒险”性(此处笔者理解有误,请看置顶回复),因为在第一个任务的频率计算后,力和位移的情况为:
Item Value Threshold Converged?
Maximum Force 0.000028 0.000450 YES
RMS Force 0.000012 0.000300 YES
Maximum Displacement 0.001345 0.001800 YES
RMS Displacement 0.000498 0.001200 YES
在第二个任务中直接读取力常数,很可能第一步仍是4个YES,导致优化直接结束。索性本例中第二个任务的第一步没有出现4个YES,优化继续下去了:
Item Value Threshold Converged?
Maximum Force 0.000028 0.000450 YES
RMS Force 0.000012 0.000300 YES
Maximum Displacement 0.116889 0.001800 NO
RMS Displacement 0.043862 0.001200 NO
实际上,在做opt+freq的任务组合时,常常出现在结构优化的最后一步是4个YES,而做频率计算后,后两个YES变成了NO的情况。这是因为在一般的结构优化中,每一步是根据近似的力常数来计算下一步的位移,而到了频率计算后,则变成了用精确的力常数来计算下一步的位移,所以后两个(或其中一个)YES经常在做了频率计算后变成了NO。一般来说,我们可以无视之,这样的结构优化结果还是可以使用的。但是如果对结果非常苛刻,则可以使用本文的方法,读取力常数,再次进行优化,或者对于不大的体系可以用上calcall或recalc=n选项。