自洽场计算的其它技巧

自洽场计算的初始猜测

自洽场计算的初始猜测轨道,对计算的收敛性有很大的影响。BDF支持多种初始猜测,如下所示:

  • Atom : 利用原子密度矩阵组合分子密度矩阵猜测,默认选项。

  • Huckel : 半经验Huckel方法猜测;

  • Hcore : 对角化单电子哈密顿猜测;

  • Readmo : 读入分子轨道做为初始猜测;

BDF默认用Atom猜测。在简洁输入模式下可以使用关键词 guess 改变BDF的初始猜测,如下所示

#! ch3cho.sh
HF/6-31G guess=Hcore unit=Bohr

geometry    # notice: unit in Bohr
C       0.1727682300       -0.0000045651       -0.8301598059
C      -2.3763311896        0.0000001634        0.5600567139
H       0.0151760290        0.0000088544       -2.9110013387
H      -2.0873396672        0.0000037621        2.5902220967
H      -3.4601725077       -1.6628370597        0.0320271859
H      -3.4601679801        1.6628382651        0.0320205364
O       2.2198078005        0.0000024315        0.2188182082
end geometry

这里,我们在第二行是用了关键词 guess=Hcore 指定使用 Hcore 猜测。SCF迭代了18次收敛。

Iter. idiis vshift  SCF Energy    DeltaE     RMSDeltaD    MaxDeltaD   Damping Times(S)
 1    0   0.000 -130.488739529 174.680929376  0.401531162  5.325668770  0.0000   0.03
 2    1   0.000 -115.595786784  14.892952744  0.407402695  5.323804678  0.0000   0.02
 3    2   0.000 -126.823748834 -11.227962049  0.115300517  1.591646800  0.0000   0.03
 4    3   0.000 -150.870636785 -24.046887951  0.011394798  0.154813426  0.0000   0.02
 5    4   0.000 -151.121829169  -0.251192384  0.004498398  0.037875784  0.0000   0.03
 6    5   0.000 -150.900123989   0.221705180  0.008483436  0.119865266  0.0000   0.02
 7    6   0.000 -151.582006133  -0.681882144  0.011892345  0.122063906  0.0000   0.02
 8    7   0.000 -152.441656890  -0.859650757  0.007907887  0.062113717  0.0000   0.03
 9    8   0.000 -152.729229838  -0.287572947  0.003318529  0.037884676  0.0000   0.02
10    2   0.000 -152.795374919  -0.066145081  0.005951772  0.054625652  0.0000   0.02
11    3   0.000 -152.839276725  -0.043901806  0.000860488  0.010210210  0.0000   0.03
12    4   0.000 -152.841131472  -0.001854746  0.000733951  0.007678730  0.0000   0.02
13    5   0.000 -152.841752921  -0.000621449  0.000348937  0.003519950  0.0000   0.02
14    6   0.000 -152.841816238  -0.000063316  0.000053288  0.000787592  0.0000   0.03
15    7   0.000 -152.841819180  -0.000002942  0.000021206  0.000157533  0.0000   0.02
16    8   0.000 -152.841819505  -0.000000325  0.000004796  0.000031694  0.0000   0.02
17    2   0.000 -152.841819522  -0.000000016  0.000000698  0.000005497  0.0000   0.03
18    3   0.000 -152.841819522  -0.000000000  0.000000236  0.000002276  0.0000   0.02
diis/vshift is closed at iter =  18
19    0   0.000 -152.8418195227 -0.000000000  0.000000078  0.000000848  0.0000   0.03

警告

这个算例分子输入坐标的单位是Bohr,必须使用关键词 unit=Bohr 指定坐标的长度单位为 Bohr

这个算例对应的高级输入为

$compass
geometry
  C 0.1727682300 -0.0000045651 -0.8301598059
  C -2.3763311896 0.0000001634 0.5600567139
  H 0.0151760290 0.0000088544 -2.9110013387
  H -2.0873396672 0.0000037621 2.5902220967
  H -3.4601725077 -1.6628370597 0.0320271859
  H -3.4601679801 1.6628382651 0.0320205364
  O 2.2198078005 0.0000024315 0.2188182082
end geometry
unit # set unit of coordinates as Bohr
  bohr
basis
  6-31g
$end

$xuanyuan
$end

$scf
rhf
guess # ask for hcore guess
  hcore
$end

备注

绝大部分情况下,Huckel和Hcore都不是最好的选择,因此除非万不得已,尽量不要使用Huckel和Hcore(尤其是后者)。

读入初始猜测轨道

BDF的SCF计算默认采用原子密度矩阵构建分子密度矩阵的方式产生初始猜测轨道。在实际计算中,用户可以读入已收敛的SCF分子轨道,做为当前SCF计算的初始猜测轨道。本算例中,我们先计算一个中性的 \(\ce{H2O}\) 分子,得到收敛轨道后,做为 \(\ce{H2O+}\) 离子的初始猜测轨道。

第一步,计算 \(\ce{H2O}\) 分子,准备输入文件,并命名为 h2o.inp 。内容如下:

#!bdf.sh
RKS/B3lyp/cc-pvdz

geometry
O
H  1  R1
H  1  R1  2 109.

R1=1.0     # OH bond length in angstrom
end geometry

执行计算后,工作目录生成可读文件 h2o.scforb ,保存了SCF计算收敛的轨道.

第二步,利用 \(\ce{H2O}\) 分子的收敛轨道做为 \(\ce{H2O+}\) 离子计算的初始猜测,准备输入文件 h2o+.inp,内容如下:

#!bdf.sh
ROKS/B3lyp/cc-pvdz guess=readmo charge=1

geometry
O
H  1  R1
H  1  R1  2 109.

R1=1.0     # OH bond length in angstrom
end geometry

%cp $BDF_WORKDIR/h2o.scforb $BDF_TMPDIR/${BDFTASK}.inporb

这里,使用了关键词 guess=readmo ,指定要读入初始猜测轨道。初始猜测轨道是用 % 引导的拷贝命令从 环境变量 BDF_WORKDIR 定义的文件夹中的h2o.scforb文件复制为 BDF_TMPDIR 中的 ${BDFTASK}.inporb 文件。 这里, BDF_WORKDIR 是执行计算任务的目录, BDF_TMPDIR 是BDF存储临时文件的目录。

与其它量子化学程序传递分子轨道

不同量子化学程序计算的分子轨道文件在原则上可以相互转化。BDF的 SCF 模块支持读入和存储 scforb 文件格式的分子轨道数据,可以通过 MOKIT(https://gitlab.com/jxzou/mokit)程序实现分子轨道文件转化,从而实现与其它量子化学程序之间传递分子轨道数据。

转化后的分子轨道文件能否正常使用,不仅依赖于原子顺序,坐标方位,以及点群对称性,还依赖于收缩函数的形式和排序,对于赝势基组还要看赝势数据是否一致。 如果有一处不一致,就会导致转化的分子轨道数据出现问题,无法达到加速收敛的目的。在基组方面,即便是同名的基组,由于原始数据的来源不同以及基组版本的差异, 计算的分子轨道数据也可能不一样,因此一定要仔细比较不同程序采用的基组是否完全一致。下面以H原子cc-pVTZ基组中的s收缩函数为例,列举了在比较基组时的几点注意事项。

****
H      1   2
S      5    3
               3.387000E+01
               5.095000E+00
               1.159000E+00
               3.258000E-01
               1.027000E-01
      6.068000E-03           0.000000E+00           0.000000E+00
      4.530800E-02           0.000000E+00           0.000000E+00
      2.028220E-01           0.000000E+00           0.000000E+00
      0.0000000E+00          1.000000E+00           0.000000E+00
      0.0000000E+00          0.000000E+00           1.000000E+00
P      2    2
()

注意

  • 收缩函数的先后顺序在两个程序中是否一致? 如果把s函数的三列收缩因子互换,分子轨道因子的排序是不一样的。

  • 每个收缩函数的收缩形式在两个程序中是否一致? 第一列收缩因子仅包含前三个s原函数,即(3s)/[1s],而在很多程序中是(4s)/[1s],那么对应的轨道因子就会有一些差别。

  • 收缩因子相位在两个程序中是否一致? 这通常出现在非收缩函数的因子 1.0 中(见第二、三列),个别程序的内置基组可能会把收缩因子 1.0 写为 -1.0,导致轨道因子差个负号。

  • 在使用赝势基组时,赝势的数据在两个程序中是否一致? 其中最典型的是 Def2系列基组的问题

为了保持基组的一致性,建议在 Compass 模块中用 ExpBas 输出基组,供其它量子化学程序使用。目前BDF支持输出 Molpro、Molcas、Gaussian、ORCA、CFour 五种基组格式。

除了以上因素之外,一些量子化学程序为了加速计算,会优先使用草稿文件夹下的临时数据文件,如果在此前的计算中恰好用了不一样的基组或者分子结构,会导致分子轨道文件非正常读取。 在BDF的计算中为了避免此类问题的发生,一般要在计算的开始把草稿文件夹清空,或者利用随机数生成新的草稿文件夹。

把小基组收敛轨道扩展为大基组初始猜测

初始猜测轨道可以由不同基组产生,同样可以加速计算收敛。这需要对初始猜测轨道文件进行扩展。 轨道扩展应该采用同组的基组,如cc-pVXZ系列、ANO-RCC系列等基组。 轨道扩展目前只支持高级输入模式。对于 \(\ce{CH3CHO}\) 分子,先用cc-pVDZ计算,然后将轨道扩展为cc-pVQZ基组计算的初始猜测轨道, 输入如下:

# First SCF calculation using small basis set cc-pvdz
$compass
geometry
C       0.1727682300       -0.0000045651       -0.8301598059
C      -2.3763311896        0.0000001634        0.5600567139
H       0.0151760290        0.0000088544       -2.9110013387
H      -2.0873396672        0.0000037621        2.5902220967
H      -3.4601725077       -1.6628370597        0.0320271859
H      -3.4601679801        1.6628382651        0.0320205364
O       2.2198078005        0.0000024315        0.2188182082
end geometry
basis
 cc-pvdz
unit # set unit of coordinates as Bohr
 Bohr
$end

$xuanyuan
$end

$scf
rhf
$end

#change chkfil name into chkfil1
%mv $BDF_WORKDIR/$BDFTASK.chkfil $BDF_WORKDIR/$BDFTASK.chkfil1

$compass
geometry
C       0.1727682300       -0.0000045651       -0.8301598059
C      -2.3763311896        0.0000001634        0.5600567139
H       0.0151760290        0.0000088544       -2.9110013387
H      -2.0873396672        0.0000037621        2.5902220967
H      -3.4601725077       -1.6628370597        0.0320271859
H      -3.4601679801        1.6628382651        0.0320205364
O       2.2198078005        0.0000024315        0.2188182082
end geometry
basis
 cc-pvqz
unit
 Bohr
$end

# change chkfil to chkfil1. notice, should use cp command since we will use
# "$BDFTASK.chkfil" in the next calculation
%cp $BDF_WORKDIR/$BDFTASK.chkfil $BDF_WORKDIR/$BDFTASK.chkfil2

# copy converged SCF orbital as input orbital of the module expandmo
%cp $BDF_WORKDIR/$BDFTASK.scforb $BDF_WORKDIR/$BDFTASK.inporb

# Expand orbital to large basis set. The output file is $BDFTASK.exporb
$expandmo
overlap
$end

$xuanyuan
$end

# use expanded orbital as initial guess orbital
%cp $BDF_WORKDIR/$BDFTASK.exporb $BDF_WORKDIR/$BDFTASK.scforb
$scf
RHF
guess
 readmo
iprtmo
 2
$end

上面的输入中,先使用 cc-pVDZ 基组执行第一个RHF计算,然后利用 expandmo 模块,将第一次 SCF 计算的收敛轨道扩展到 cc-pVQZ 基组, 最后利用 guess=readmo 做为SCF要读入的初始猜测轨道。

expandmo模块的输出为,

|******************************************************************************|

    Start running module expandmo
    Current time   2021-11-29  22:20:50

|******************************************************************************|
 $expandmo
 overlap
 $end
 /Users/bsuo/check/bdf/bdfpro/ch3cho_exporb.chkfil1
 /Users/bsuo/check/bdf/bdfpro/ch3cho_exporb.chkfil2
 /Users/bsuo/check/bdf/bdfpro/ch3cho_exporb.inporb
  Expanding MO from small to large basis set or revise ...

 1 Small basis sets

 Number of  basis functions (NBF):      62
 Maxium NBF of shell :        6

 Number of basis functions of small basis sets:       62

 2 Large basis sets

 Number of  basis functions (NBF):     285
 Maxium NBF of shell :       15

  Overlap expanding :                     1
 Read guess orb
 Read orbital title:  TITLE - SCF Canonical Orbital
nsbas_small  62
nsbas_large 285
ipsmall   1
iplarge   1
  Overlap of dual basis ...
  Overlap of large basis ...
 Write expanded MO to scratch file ...
|******************************************************************************|

    Total cpu     time:          0.42  S
    Total system  time:          0.02  S
    Total wall    time:          0.47  S

    Current time   2021-11-29  22:20:51
    End running module expandmo
|******************************************************************************|

可以看出,小基组有62个轨道,大基组有285个轨道,expandmo读入了SCF收敛的正则轨道,扩展到大基组并写入临时文件。

第二次SCF计算的输出为,

  /Users/bsuo/check/bdf/bdfpro/ch3cho_exporb.scforb
  Read guess orb:  nden=1  nreps= 1  norb=  285  lenmo=  81225
  Read orbital title:  TITLE - orthognal Expand CMO
  Orbitals initialization is completed.

  ........
Iter. idiis vshift  SCF Energy    DeltaE     RMSDeltaD    MaxDeltaD   Damping Times(S)
 1    0   0.000 -152.952976892 122.547522034  0.002218985  0.246735859  0.0000  16.30
 2    1   0.000 -152.983462881  -0.030485988  0.000367245  0.026196100  0.0000  16.83
 3    2   0.000 -152.983976045  -0.000513164  0.000086429  0.006856831  0.0000  17.18
 4    3   0.000 -152.984012062  -0.000036016  0.000016763  0.001472939  0.0000  17.02
 5    4   0.000 -152.984019728  -0.000007666  0.000010400  0.001012788  0.0000  17.42
 6    5   0.000 -152.984021773  -0.000002045  0.000003396  0.000328178  0.0000  17.28
 7    6   0.000 -152.984022197  -0.000000423  0.000001082  0.000075914  0.0000  17.40
 8    7   0.000 -152.984022242  -0.000000044  0.000000154  0.000008645  0.0000  17.28
 9    8   0.000 -152.984022243  -0.000000001  0.000000066  0.000005087  0.0000  19.38
diis/vshift is closed at iter =   9
10    0   0.000 -152.984022243  -0.000000000  0.000000007  0.000000584  0.0000  18.95

    Label              CPU Time        SYS Time        Wall Time
   SCF iteration time:       517.800 S        0.733 S      175.617 S

收敛到具有特定自旋布居的SCF波函数

对于某些电子结构复杂的体系,尤其是存在反铁磁耦合的过渡金属配合物体系,默认的初猜常常无法收敛到正确的SCF解。典型例子如含两个Fe(III)的铁硫簇模型体系[Fe2S2(SR)4]2-,基态为单重态,可能收敛到(1)闭壳层单重态;(2)两个铁均为低自旋(S=1/2),自旋方向相反的反铁磁耦合态;(3)两个铁均为低自旋(S=3/2),自旋方向相反的反铁磁耦合态;(4)两个铁均为高自旋(S=5/2),自旋方向相反的反铁磁耦合态。直接做RKS计算会收敛到(1);以三重态为初猜计算单重态,虽往往能打破自旋对称性,但一般会收敛到(2);而体系的真正基态为(4)。为了收敛到(4),可以采用 FLMO方法 ,先将体系分为两个[Fe(SR)2]+片段和两个S2-片段,分别计算得到片段波函数后再组装得到总体系波函数,但这需要用户进行较多操作。

针对FLMO方法使用不便或难以分片的情形,自2025年8月起,BDF还提供一种不借助分片方法即可收敛到指定SCF解的方法,只需直接在 $scf 块中用 restrainspin 关键词指定哪些原子需要占据几个单电子即可。该方法正确收敛的概率并非100%,部分原因在于用户指定的单电子占据情况可能并不对应一个稳定的SCF解,但出于同样的原因,FLMO及其他方法的收敛概率也无法达到100%。例如以下算例计算含一个Fe(III)、一个Fe(II)的簇[Fe2S2(SH)4]3-,且要求两个Fe原子都是高自旋、彼此反铁磁耦合:

$compass
title
 dinuclear Fe complex
basis-block
 3-21G # for real calculations, use at least 6-31G(d)
 Fe = LANL2DZ # for real calculations, def2-TZVP is recommended
end basis
geometry
Fe     14.5260939424     12.5377828892      2.1988080086
 S     16.1083103440     11.1748725273      1.5865791609
 S     12.9003852428     11.7308390465      1.1881934067
Fe     14.4298509469     10.6877447671      0.4025405615
 S     14.3280639036     13.3338865932      4.3110023010
 S     15.5099134429     14.4139310766      1.7666541957
 S     13.7633959005      8.7510441770      1.5795853408
 S     14.1408273025     10.0269442277     -1.6706111612
 H     13.2250185390      9.1777476902     -2.1406096610
 H     16.7916110121     14.2772746738      2.1363356902
 H     14.2101464458     12.1727249125      4.9478734403
 H     12.9707766368      8.4446622752      0.5668718853
end geometry
norotate
mpec+cosx
$end

$xuanyuan
$end

$scf
uks
dft
 TPSSh
charge
 -3
spinmulti
 2
restrainspin
 2 # restrain the spin populations of two atoms
 1 +5.0 # the 1st atom (Fe) should have 5 alpha spins, i.e. high-spin Fe(III)
 4 -4.0 # the 4th atom (Fe) should have 4 beta spins, i.e. high-spin Fe(II)
maxiter
 300
vshift
 0.5
damp
 0.7
grid
 fine
solvent
 water
$end

计算收敛到能量为-2627.17018238 Hartree的一个态,且Mulliken自旋布居符合要求:

[Mulliken Population Analysis]
 Atomic charges and Spin densities :
    1Fe     -0.1765    3.2812
    2S      -0.4581    0.2066
    3S      -0.3378    0.1913
    4Fe     -0.2139   -2.9579
    5S      -0.4816    0.1892
    6S      -0.4496    0.2777
    7S      -0.5426   -0.1005
    8S      -0.5923   -0.1063
    9H       0.0595    0.0085
   10H       0.0888    0.0080
   11H       0.0546   -0.0023
   12H       0.0497    0.0044
    Sum:    -3.0000    1.0000

这里两个Fe的自旋布居(3左右)比预期值小是正常现象,因为有一部分自旋密度离域到了附近的原子上,如需进一步确认,可作前线轨道图辅助判断。若去掉 restrainspin 及后面的3行,程序将收敛到能量较高的一个态(-2627.12435953 Hartree),其中两个Fe都是低自旋:

[Mulliken Population Analysis]
 Atomic charges and Spin densities :
    1Fe     -0.3877   -0.2409
    2S      -0.4207    0.0492
    3S      -0.2650   -0.0768
    4Fe     -0.2587    1.2475
    5S      -0.4607   -0.0040
    6S      -0.4424    0.0551
    7S      -0.4617   -0.0417
    8S      -0.5385    0.0338
    9H       0.0614    0.0005
   10H       0.0700   -0.0251
   11H       0.0470    0.0029
   12H       0.0570   -0.0006
    Sum:    -3.0000    1.0000

由此可以看出, restrainspin 对于反铁磁耦合体系收敛到正确的解非常有帮助。

分子轨道最大占据数(mom)方法计算激发态

mom(maximum occupation method)是一种ΔSCF方法,可用于计算激发态。注意该方法的缩写为全小写字母,以和另一种ΔSCF方法——MOM(maximum overlap method)方法区分。

#----------------------------------------------------------------------
#
# mom method: J. Liu, Y. Zhang, and W. Liu, J. Chem. Theory Comput. 10, 2436 (2014).
#
# gs  = -169.86584128
# ab  = -169.62226127
# T   = -169.62483480
# w(S)= 6.69eV
#----------------------------------------------------------------------
$COMPASS
Title
 mom
Basis
 6-311++GPP
Geometry
 C       0.000000    0.418626    0.000000
 H      -0.460595    1.426053    0.000000
 O       1.196516    0.242075    0.000000
 N      -0.936579   -0.568753    0.000000
 H      -0.634414   -1.530889    0.000000
 H      -1.921071   -0.362247    0.000000
End geometry
Check
$END

$XUANYUAN
$END

$SCF
UKS
DFT
B3LYP
alpha
  10 2
beta
  10 2
$END

%cp ${BDFTASK}.scforb $BDF_TMPDIR/${BDFTASK}.inporb

# delta scf with mom
$SCF
UKS
DFT
B3LYP
guess
 readmo
alpha
 10 2
beta
 10 2
ifpair
hpalpha
 1
 10 0
 11 0
iaufbau
 2
$END

# pure delta scf for triplet
$SCF
UKS
DFT
B3LYP
alpha
  11 2
beta
  9 2
$END

这个算例执行了三次SCF计算,

  • 第一次SCF,利用UKS方法计算甲酰胺分子的基态。输入利用alpha与beta两个关键词,分别指定了alpha和beta轨道的占据情况。甲酰胺分子基态是单重态S0,这里指定的alpha和beta占据情况相同。 10 2 指定不可约表示A'与A"分别有10个和2个占据轨道。SCF模块将根据构造原理,按照轨道能量由低到高填充电子到轨道上。

  • 第二次SCF,利用UKS与mom方法计算甲酰胺分子的S1态。这里的关键点有:1 利用guess=readmo指定读入上一步UKS的收敛轨道;2 利用alpha、beta关键词设置了每个对称性轨道的占据数;3 设置了变量ifpair,需要和hpalpha,hpbeta联用,用于指定空穴-粒子(hole-particle - HP)轨道对的电子激发情况;4 设置了hpalpha变量,指定激发的HP轨道对。数字1表示激发一对HP轨道,下面的两行指定轨道激发情况,第一列表示在第一个不可约表示中把第10个alpha轨道的电子激发到第11个alpha轨道;第二列元素都为零,表示第二个不可约表示的轨道不做激发; 5 iaufbau变量设置为2,指定要进行mom计算。

  • 第三次SCF,利用UKS方法计算甲酰胺分子的T1态。输入中,我们利用alpha和beta关键词指定轨道占据情况,其中alpha轨道的占据数为 11 2 ,表示对称性为A'和A"的alpha轨道上分别有11和2个电子占据, beta轨道的占据情况为 9 2 。 因为所要求解的态是给定的轨道占据数下能量最低的态,因此无需指定iaufbau。

这里,第一次SCF计算收敛结果为,

 Superposition of atomic densities as initial guess.
 skipaocheck T F
 Solve HC=EC in pflmo space. F       12       75
 Initial guess energy =   -169.2529540680

 [scf_cycle_init_ecdenpot]
Meomory for coulpotential         0.00  G

 Start SCF iteration......

Iter. idiis vshift  SCF Energy    DeltaE     RMSDeltaD    MaxDeltaD   Damping Times(S)
 1    0   0.000 -169.411739263  -0.158785195  0.005700928  0.163822560  0.0000   0.20
Turn on DFT calculation ...
 2    1   0.000 -169.743175119  -0.331435856  0.008905349  0.340815886  0.0000   0.42
 3    2   0.000 -169.232333660   0.510841459  0.006895796  0.296788710  0.0000   0.43
 4    3   0.000 -169.863405142  -0.631071482  0.000364999  0.015732911  0.0000   0.43
 5    4   0.000 -169.863345847   0.000059295  0.000209771  0.009205878  0.0000   0.42
 6    5   0.000 -169.865811301  -0.002465454  0.000027325  0.000606909  0.0000   0.43
 7    6   0.000 -169.865831953  -0.000020651  0.000008039  0.000357726  0.0000   0.43
 8    7   0.000 -169.865833199  -0.000001246  0.000003927  0.000114311  0.0000   0.42
 9    8   0.000 -169.865833401  -0.000000201  0.000000182  0.000004399  0.0000   0.43
diis/vshift is closed at iter =   9
10    0   0.000 -169.865833402  -0.000000000  0.000000139  0.000003885  0.0000   0.43

  Label              CPU Time        SYS Time        Wall Time
 SCF iteration time:         8.650 S        0.700 S        4.050 S

 Final DeltaE =  -4.4343551053316332E-010
 Final DeltaD =   1.3872600382452641E-007   5.0000000000000002E-005

 Final scf result
   E_tot =              -169.86583340
   E_ele =              -241.07729109
   E_nn  =                71.21145769
   E_1e  =              -371.80490197
   E_ne  =              -541.14538673
   E_kin =               169.34048477
   E_ee  =               148.48285541
   E_xc  =               -17.75524454
  Virial Theorem      2.003102

可以看出,第一次SCF计算使用了atom猜测,计算得到S0的能量为 -169.8658334023 a.u. 。第二次SCF计算读入了第一次SCF的收敛轨道, 并使用mom方法做SCF计算,输出文件先提示读入了分子轨道,并给出占据情况,

   [Final occupation pattern: ]

Irreps:        A'      A''

detailed occupation for iden/irep:      1   1
   1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00
detailed occupation for iden/irep:      1   2
   1.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00
Alpha      10.00    2.00

这里, A' 不可约表示的第10个alpha轨道是占据轨道,第11个轨道是空轨道。第二次SCF计算读入了第一次SCF的收敛轨道,并使用mom方法做SCF计算,输入中要求将 A' 表示的第10个轨道的电子激发到第11个轨道上。输出文件先提示读入了分子轨道,并给出占据情况,

Read initial orbitals from user specified file.

/tmp/20117/mom_formamide.inporb
Read guess orb:  nden=2  nreps= 2  norb=   87  lenmo=   4797
Read orbital title:  TITLE - SCF Canonical Orbital

Initial occupation pattern: iden=1  irep= 1  norb(irep)=   66
   1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.00
   1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00


Initial occupation pattern: iden=1  irep= 2  norb(irep)=   21
   1.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00


Initial occupation pattern: iden=2  irep= 1  norb(irep)=   66
   1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00


Initial occupation pattern: iden=2  irep= 2  norb(irep)=   21
   1.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
   0.00

这里,iden=1为alpha轨道,irep=1指第一个不可约表示,总共有norb=66个轨道,其中,第10个轨道的占据数为0.00,第11个轨道占据数为1.00。经14次SCF迭代,收敛的S1态能量为 -169.6222628003 a.u.,如下所示:

Iter. idiis vshift  SCF Energy    DeltaE     RMSDeltaD    MaxDeltaD   Damping Times(S)
 1    0   0.000 -169.505632070 125.031578610  0.020428031  1.463174456  0.0000   0.45
 2    1   0.000 -169.034645773   0.470986296  0.036913522  1.562284831  0.0000   0.43
 3    2   0.000 -165.750862892   3.283782881  0.032162782  1.516480990  0.0000   0.43
 4    3   0.000 -169.560678610  -3.809815718  0.008588866  0.807859419  0.0000   0.43
 5    4   0.000 -169.596211021  -0.035532411  0.003887621  0.367391029  0.0000   0.42
 6    5   0.000 -169.620128518  -0.023917496  0.001826050  0.172456003  0.0000   0.43
 7    6   0.000 -169.621976725  -0.001848206  0.000486763  0.044630527  0.0000   0.43
 8    7   0.000 -169.622245116  -0.000268391  0.000113718  0.004980035  0.0000   0.43
 9    8   0.000 -169.622261269  -0.000016153  0.000112261  0.009715905  0.0000   0.42
10    2   0.000 -169.622262553  -0.000001284  0.000043585  0.004092668  0.0000   0.42
11    3   0.000 -169.622262723  -0.000000169  0.000031601  0.002792075  0.0000   0.42
12    4   0.000 -169.622262790  -0.000000067  0.000010125  0.000848297  0.0000   0.43
13    5   0.000 -169.622262798  -0.000000007  0.000003300  0.000273339  0.0000   0.43
 diis/vshift is closed at iter =  13
14    0   0.000 -169.622262800  -0.000000002  0.000001150  0.000079378  0.0000   0.42

  Label              CPU Time        SYS Time        Wall Time
 SCF iteration time:        13.267 S        0.983 S        6.000 S

 Final DeltaE =  -1.8403909507469507E-009
 Final DeltaD =   1.1501625138328933E-006   5.0000000000000002E-005

 Final scf result
   E_tot =              -169.62226280
   E_ele =              -240.83372049
   E_nn  =                71.21145769
   E_1e  =              -368.54021347
   E_ne  =              -537.75897296
   E_kin =               169.21875949
   E_ee  =               145.28871749
   E_xc  =               -17.58222451
  Virial Theorem      2.002385


 [Final occupation pattern: ]

 Irreps:        A'      A''

 detailed occupation for iden/irep:      1   1
    1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.00
    1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
    0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
    0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
    0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
    0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
    0.00 0.00 0.00 0.00 0.00 0.00

SCF收敛后再一次打印轨道占据情况,可以看到 alpha 轨道中 A' 不可约表示的第10个轨道没有电子占据,第11个轨道有一个电子占据。

第三个SCF计算给出了 T1 态能量,为 -169.6248370697 a.u.,输出如下:

Iter. idiis vshift  SCF Energy    DeltaE     RMSDeltaD    MaxDeltaD   Damping Times(S)
 1    0   0.000 -169.411739263  -0.158785195  0.083821477  9.141182225  0.0000   0.17
 Turn on DFT calculation ...
 2    1   0.000 -169.480549474  -0.068810211  0.066700318  6.978728919  0.0000   0.40
 3    2   0.000 -169.277735673   0.202813801  0.014778190  0.648183923  0.0000   0.42
 4    3   0.000 -169.613991196  -0.336255522  0.005923909  0.621843348  0.0000   0.42
 5    4   0.000 -169.620096778  -0.006105582  0.001967168  0.164506160  0.0000   0.40
 6    5   0.000 -169.623636999  -0.003540220  0.002722812  0.246425639  0.0000   0.42
 7    6   0.000 -169.624704514  -0.001067515  0.001064536  0.098138798  0.0000   0.42
 8    7   0.000 -169.624814882  -0.000110368  0.000525436  0.046392861  0.0000   0.42
 9    8   0.000 -169.624834520  -0.000019637  0.000179234  0.012966641  0.0000   0.42
10    2   0.000 -169.624836694  -0.000002174  0.000063823  0.004902276  0.0000   0.42
11    3   0.000 -169.624836922  -0.000000227  0.000017831  0.001440089  0.0000   0.43
12    4   0.000 -169.624837025  -0.000000103  0.000034243  0.002618897  0.0000   0.42
13    5   0.000 -169.624837065  -0.000000039  0.000006158  0.000466001  0.0000   0.40
14    6   0.000 -169.624837068  -0.000000003  0.000003615  0.000354229  0.0000   0.42
diis/vshift is closed at iter =  14
15    0   0.000 -169.624837069  -0.000000001  0.000000966  0.000070404  0.0000   0.42

 Label              CPU Time        SYS Time        Wall Time
SCF iteration time:        13.150 S        0.950 S        5.967 S

Final DeltaE =  -1.1375220765330596E-009
Final DeltaD =   9.6591808698539483E-007   5.0000000000000002E-005

Final scf result
  E_tot =              -169.62483707
  E_ele =              -240.83629476
  E_nn  =                71.21145769
  E_1e  =              -368.57834907
  E_ne  =              -537.80483706
  E_kin =               169.22648799
  E_ee  =               145.32683246
  E_xc  =               -17.58477815
 Virial Theorem      2.002354

备注

对于某些体系的mom计算,BDF默认开启的SMH收敛算法可能反倒会阻碍收敛,此时可尝试在$scf块中加入NoSMH关键字,有一定概率会使得SCF收敛。如仍不奏效,可参照下节解决SCF不收敛问题的方法,来解决mom计算的收敛问题。

处理自洽场计算的不收敛问题

当SCF计算完成后,用户务必检查SCF是否收敛,只有在收敛的前提下才可以使用SCF计算的结果(能量,布居分析,轨道能等)以及进行后续的计算。注意SCF是否收敛不能仅从输出文件末尾有没有报错来判断,因为即便SCF不收敛,程序也不会立刻退出,而只是在SCF迭代的输出之后、SCF能量的输出之前,提示:

Warning !!! Total energy not converged!

而即便在这种情况下,程序仍然会在该信息之后打印能量、轨道信息、布居分析结果等,其中SCF能量后面会有 (NOT CONVERGED) 字样。虽然这些结果不能作为正式计算结果使用,但它们对于分析SCF不收敛的原因有一定帮助。

导致SCF不收敛的常见原因包括:

  1. HOMO-LUMO能隙过小,导致前线轨道的占据情况反复变化。例如两个轨道 \(\psi_1\)\(\psi_2\) ,在第N次SCF迭代时 \(\psi_1\) 为占据轨道, \(\psi_2\) 为空轨道,然而以这样的轨道占据情况为基础构建Fock矩阵并对角化后,得到的第N+1次SCF迭代的轨道,却是 \(\psi_1\) 的轨道能较 \(\psi_2\) 更高,因此电子从 \(\psi_1\) 轨道转移到 \(\psi_2\) 轨道。但这样一来,第N+1次SCF迭代的Fock矩阵相比第N次SCF迭代就会发生很大变化,导致在第N+2次SCF迭代时 \(\psi_1\) 的轨道能较 \(\psi_2\) 更低,于是轨道占据数又回到了第N次SCF迭代的情形,因而SCF迭代的轨道占据数总是在变化,始终不收敛。这种情况的典型表现为SCF能量交替在两个能量之间振荡(或在一定范围内无规律振荡),振荡幅度在 \(10^{-4} \sim 1\) Hartree左右,且SCF结束后打印的轨道占据数与预期不符。

  2. HOMO-LUMO能隙较小,虽然各步迭代的轨道占据数没有变化,但轨道形状反复变化,导致SCF振荡不收敛。这种情况的典型表现与前一条类似,但振荡幅度一般稍小些,且SCF结束后打印的轨道占据数与预期定性相符。

  3. 数值积分格点过小或者双电子积分精度过低,导致SCF因数值误差而小幅度振荡不收敛。这种情况的典型表现为SCF能量以 \(10^{-4}\) Hartree以下的幅度无规律振荡,且SCF结束后打印的轨道占据数与预期定性相符。

  4. 基组接近线性相关,或因为格点太小导致基组在格点上的投影接近线性相关。这种情况的典型表现为SCF能量以1 Hartree以上的幅度变化(不一定是在振荡,也可能是单调或者基本单调的变化),SCF能量远低于预期值,且SCF结束后打印的轨道占据数完全不符合物理实际。当SCF能量较预期值低得非常多时,SCF能量甚至可能不显示为数字,而是显示为一串星号。

以下是各类SCF不收敛问题的常见解决方法(一定程度上也适用于BDF以外的软件):

  1. 增加能级移动vshift,适用于第1类和第2类情况,方法为在输入文件的$scf模块里加入:

vshift
 0.2

如果仍然观察到明显的振荡,则逐渐增加vshift,直到收敛为止。vshift会倾向于让SCF的收敛变得单调,但是vshift设得太大会增加迭代收敛所用的次数,因此在增加vshift的时候可以适当增加maxiter。当vshift增加到1.0仍然无法收敛时,应该考虑其他方法。

  1. 增加密度矩阵阻尼damp,适用于第2类情况(对第1类情况也有一点效果),方法为在输入文件的$scf模块里加入:

damp
 0.7

注意damp可以和vshift联用,两者的效果在一定程度上是相互促进的。如果阻尼设为0.7仍然观察到明显的振荡,则在保证阻尼小于1的情况下增大阻尼,例如接下来可以尝试0.9、0.95等。与vshift类似,damp也是倾向于改善SCF收敛的单调性,但damp太大会导致收敛变慢,因此可以增加maxiter。当damp增加到0.99仍然无法收敛时,应该考虑其他方法。

  1. 关闭DIIS,适用于第1类和第2类情况,且增加vshift和damp也无法收敛时。DIIS在大多数情况下是会加速SCF收敛的,但当HOMO-LUMO能隙特别小时有可能反倒会减慢甚至阻止收敛,后一种情况下可以在$scf模块里添加NoDIIS关键词关掉DIIS,增加maxiter,并视收敛情况设定vshift和damp。

  2. 关闭SMH,适用于第1类和第2类情况,且前3种方式都不奏效时,方法是在$scf模块里添加NoSMH关键词,增加maxiter,并视收敛情况设定vshift和damp。目前看来,起码对于基态计算而言,用SMH不收敛、不用SMH反倒能收敛的情形极其少,但是因为SMH是一种很新的SCF收敛方法 [66] ,不排除极个别情况下SMH会对收敛有负面影响,因此关闭SMH可以作为一种备选方案。

  3. 改用FLMO或iOI方法,适用于第1类和第2类情况,分子较大(如大于50个原子),且怀疑SCF不收敛是因为原子初始猜测精度太低或者定性错误所导致时。方法请参见 FLMO及iOI方法相关章节

  4. 先计算一个类似的、较容易收敛的体系,再以该体系的波函数为初猜来收敛原体系,适用于第1类和第2类情况。比如一个中性的二重态过渡金属配合物的SCF计算不收敛,可以计算其闭壳层的一价阳离子,收敛后以一价阳离子的轨道作为初猜来进行中性分子的SCF计算(但注意因为BDF尚不支持读取RHF/RKS波函数作为UHF/UKS计算的初猜,因此此处闭壳层的一价阳离子应当用UHF/UKS计算)。极端情况下甚至可以先计算高价阳离子,然后添加少量(如2个)电子重新收敛SCF,再添加少量电子,如此直至得到原来的中性体系的波函数。另一种常用的手段为先在小基组下进行SCF计算,收敛后利用 expandmo模块 将小基组的SCF轨道投影到原基组上,再在原基组下进行SCF迭代直至收敛。

  5. 增大格点,适用于第3类情况,有时对第4类情况也有效。方法是用grid关键词,如:

grid
 fine

注意:(1)对于meta-GGA泛函,默认的格点已经是fine了,因此此时应当将格点设为ultra fine;(2)增大格点会增加每一步SCF迭代的耗时;(3)增大格点会使得收敛的能量和其他没有改变grid的计算不可比,因此如果要将这个计算和以前做过的计算进行比较,或者将这个计算得到的能量/自由能与其他计算的结果作差等等,则必须把已经做过的所有相关计算用和本输入文件相同的格点重新计算,即便已经做过的那些计算不加大格点也能收敛,也需要这样做。加大格点后若结果没有任何改善,则应该尝试其他方法;如果结果有改善但还是不收敛,可以进一步尝试将fine改为ultra fine;如果仍然不能收敛,应当考虑下面的方法。

  1. 关闭incremental Fock,适用于第3类情况。方法是在SCF模块里添加NoDeltaP。这会增加每一步SCF迭代的耗时,但结果和不加NoDeltaP的结果是可比的。

  2. 将双电子积分的阈值设严,适用于第3类情况,有时对第4类情况也有效。方法是在SCF模块里添加:

optscreen
 1

该方法和增大格点一样,也会增大每一步SCF迭代的耗时,且也会导致计算结果和不加optscreen的计算结果不可比。该方法仅适用于不开启MPEC或MPEC+COSX的计算。

  1. 将判断基组线性相关性的阈值设松,适用于第4类情况。方法是在$scf模块里添加:

checklin
tollin
 1.d-6

该方法会导致计算结果和不加这些关键词的计算结果不可比。tollin不建议设得比1.d-5更大,否则会引入较大误差,如果tollin设为1.d-5仍然出现第4类不收敛情况,则应考虑以上所述的增大格点、改变双电子积分阈值等方法。

注意在以上各方法中,如果某种方法虽不能使SCF收敛,但让SCF收敛情况较以前更好了,则尝试下一个方法时应当用

guess
 readmo

读取上一种方法的最后一步SCF迭代的轨道作为初猜。但如果前一种方法反倒导致SCF收敛变差了,则尝试下一个方法时应当重新从原子猜测开始,或者挑选之前尝试过的其他方法的最后一步迭代的轨道作为初猜(当然这要求用户提前把每种SCF收敛方法得到的轨道都进行备份)。

自洽场计算的加速算法

BDF的一个重要特色是利用 MPEC+COSX 方法加速SCF、TDDFT的能量及梯度计算。设置MPEC+COSX计算,输入如下:

#! amylose2.sh
HF/cc-pvdz  MPEC+COSX

Geometry
H      -5.27726610038004     0.15767995434597     1.36892178079618
H      -3.89542800401751    -2.74423996083456    -2.30130324998720
H      -3.40930212959730     3.04543096108345     1.73325487719318
O      -4.25161610042910    -0.18429704053319     1.49882079466485
H      -4.12153806480025     0.39113300040060    -0.47267019103680
O      -3.93883902709049    -2.16385597983528    -1.37984323910654
H      -3.65755506365314    -2.55190701717719     0.56784675873394
H      -2.66688104102718    -3.13999999152083    -0.32869523309397
O      -3.68737510690803     2.57255697808269     0.79063986197194
H      -2.16845111442446     1.40439897322928     1.59675986910159
H      -0.80004208156425     3.67692503357694    -0.87083105709857
C      -3.47036908085237     0.21757398797107     0.38361581865084
C      -3.08081604941874    -2.23618399620817    -0.25179522317288
H      -1.85215308213129    -1.05270701067006     0.92020982572454
C      -2.73634509645279     1.50748698767418     0.67208385967460
O      -0.95388209186676     2.93603601652216    -0.08659407523165
H      -2.34176605974133     2.08883703173396    -1.35500112054343
C      -2.46637306624908    -0.89337899823852     0.07760781649778
C      -1.77582007601201     1.83730601785282    -0.45887211416401
O      -1.70216504605578    -0.48600696920536    -1.07005315975028
H      -0.26347504436884     0.90841605388912    -1.67304510231922
C      -0.87599906000257     0.65569503172715    -0.80788211986139
H       1.05124197574425    -4.08129295376550    -0.80486617677089
H       1.91283792081157     2.93924205088598    -0.71300301703422
O       0.07007992244287     0.29718501862843     0.19143889205868
H       1.28488995808993    -0.48228594245462    -1.27588009910221
O       0.83243796215244    -3.05225096122844    -0.51820416035526
H       0.03099092283770    -2.15700599981123     1.08682384153403
H       0.99725792474852    -3.21082099855794     1.38542783977374
O       1.92550793896406     1.99389906198042    -1.25576903593383
H       2.32288890226196     1.52348902475463     0.72949896259198
H       5.42304993860699     1.71940008598879    -1.13583497057179
C       1.35508593454345    -0.11004196264200    -0.25348109013556
C       0.98581793175676    -2.43946398581436     0.75228585517262
H       1.91238990103301    -0.83125899736406     1.66788890655085
C       2.32240292575108     1.05122704465611    -0.25278704698785
O       4.65571492366175     1.63248206459704    -0.36643098789343
H       3.77658595927138     0.23304608296485    -1.60079803407907
C       1.86060292384221    -1.20698497780059     0.68314589788694
C       3.72997793572998     0.57134806164321    -0.56599702816882
O       3.14827793673614    -1.62888795836893     0.20457391544942
H       5.12279093584136    -0.96659193933436     0.00181296891020
C       4.14403492674986    -0.60389595307832     0.31494395641232
O       4.31314989648861    -0.29843197973243     1.69336596603165
H       3.37540288537848     0.07856300492440     2.10071295465512
End geometry

如果在高级输入模式下,只需在COMPASS模块输入中加入关键词 MPEC+COSX,如:

$compass
Geometry
H      -5.27726610038004     0.15767995434597     1.36892178079618
H      -3.89542800401751    -2.74423996083456    -2.30130324998720
H      -3.40930212959730     3.04543096108345     1.73325487719318
O      -4.25161610042910    -0.18429704053319     1.49882079466485
H      -4.12153806480025     0.39113300040060    -0.47267019103680
O      -3.93883902709049    -2.16385597983528    -1.37984323910654
H      -3.65755506365314    -2.55190701717719     0.56784675873394
H      -2.66688104102718    -3.13999999152083    -0.32869523309397
O      -3.68737510690803     2.57255697808269     0.79063986197194
H      -2.16845111442446     1.40439897322928     1.59675986910159
H      -0.80004208156425     3.67692503357694    -0.87083105709857
C      -3.47036908085237     0.21757398797107     0.38361581865084
C      -3.08081604941874    -2.23618399620817    -0.25179522317288
H      -1.85215308213129    -1.05270701067006     0.92020982572454
C      -2.73634509645279     1.50748698767418     0.67208385967460
O      -0.95388209186676     2.93603601652216    -0.08659407523165
H      -2.34176605974133     2.08883703173396    -1.35500112054343
C      -2.46637306624908    -0.89337899823852     0.07760781649778
C      -1.77582007601201     1.83730601785282    -0.45887211416401
O      -1.70216504605578    -0.48600696920536    -1.07005315975028
H      -0.26347504436884     0.90841605388912    -1.67304510231922
C      -0.87599906000257     0.65569503172715    -0.80788211986139
H       1.05124197574425    -4.08129295376550    -0.80486617677089
H       1.91283792081157     2.93924205088598    -0.71300301703422
O       0.07007992244287     0.29718501862843     0.19143889205868
H       1.28488995808993    -0.48228594245462    -1.27588009910221
O       0.83243796215244    -3.05225096122844    -0.51820416035526
H       0.03099092283770    -2.15700599981123     1.08682384153403
H       0.99725792474852    -3.21082099855794     1.38542783977374
O       1.92550793896406     1.99389906198042    -1.25576903593383
H       2.32288890226196     1.52348902475463     0.72949896259198
H       5.42304993860699     1.71940008598879    -1.13583497057179
C       1.35508593454345    -0.11004196264200    -0.25348109013556
C       0.98581793175676    -2.43946398581436     0.75228585517262
H       1.91238990103301    -0.83125899736406     1.66788890655085
C       2.32240292575108     1.05122704465611    -0.25278704698785
O       4.65571492366175     1.63248206459704    -0.36643098789343
H       3.77658595927138     0.23304608296485    -1.60079803407907
C       1.86060292384221    -1.20698497780059     0.68314589788694
C       3.72997793572998     0.57134806164321    -0.56599702816882
O       3.14827793673614    -1.62888795836893     0.20457391544942
H       5.12279093584136    -0.96659193933436     0.00181296891020
C       4.14403492674986    -0.60389595307832     0.31494395641232
O       4.31314989648861    -0.29843197973243     1.69336596603165
H       3.37540288537848     0.07856300492440     2.10071295465512
End geometry
Basis
  cc-pvdz
MPEC+COSX # ask for the MPEC+COSX method
$end

SCF 模块会输出会有关 MPEC+COSX 是否都被设置为 True 的提示:

--- PRINT: Information about SCF Calculation ---
ICTRL_FRAGSCF=  0
IPRTMO=  1
MAXITER=  100
THRENE= 0.10E-07 THRDEN= 0.50E-05
DAMP= 0.00 VSHIFT= 0.00
IFDIIS= T
THRDIIS= 0.10E+01
MINDIIS=   2 MAXDIIS=   8
iCHECK=  0
iAUFBAU=  1
INIGUESS=  0
IfMPEC= T
IfCOSX= T

这里, IfMPEC= T , 且 IfCOSX= T 说明 MPEC+COSX 方法被用于计算。SCF迭代过程如下:

 [scf_cycle_init_ecdenpot]
Meomory for coulpotential         0.02  G

 Start SCF iteration......


Iter.   idiis  vshift       SCF Energy            DeltaE          RMSDeltaD          MaxDeltaD      Damping    Times(S)
   1      0    0.000   -1299.6435521238     -23.7693069405       0.0062252375       0.2842668435    0.0000      2.69
   2      1    0.000   -1290.1030630508       9.5404890730       0.0025508000       0.1065204344    0.0000      1.65
   3      2    0.000   -1290.2258798561      -0.1228168053       0.0014087449       0.0742227520    0.0000      1.67
   4      3    0.000   -1290.4879683983      -0.2620885422       0.0002338141       0.0153879051    0.0000      1.64
   5      4    0.000   -1290.4955210658      -0.0075526675       0.0000713807       0.0049309441    0.0000      1.57
   6      5    0.000   -1290.4966349620      -0.0011138962       0.0000156009       0.0010663736    0.0000      1.51
   7      6    0.000   -1290.4966797420      -0.0000447800       0.0000043032       0.0002765334    0.0000      1.44
   8      7    0.000   -1290.4966810419      -0.0000012999       0.0000014324       0.0000978302    0.0000      1.37
   9      8    0.000   -1290.4966794202       0.0000016217       0.0000003030       0.0000173603    0.0000      1.40
  10      2    0.000   -1290.4966902283      -0.0000108081       0.0000000659       0.0000034730    0.0000      1.11
 diis/vshift is closed at iter =  10
  11      0    0.000   -1290.5003691464      -0.0036789181       0.0000225953       0.0009032949    0.0000      5.85

  Label              CPU Time        SYS Time        Wall Time
 SCF iteration time:       179.100 S        1.110 S       22.630 S

 Final DeltaE = -3.678918126752251E-003
 Final DeltaD =  2.259533940614071E-005  5.000000000000000E-005

 Final scf result
   E_tot =             -1290.50036915
   E_ele =             -3626.68312754
   E_nn  =              2336.18275840
   E_1e  =             -6428.96436179
   E_ne  =             -7717.90756825
   E_kin =              1288.94320647
   E_ee  =              2802.28123424
   E_xc  =                 0.00000000
  Virial Theorem      2.001208

在CPU为i9-9900K的台式机上,8个OpenMP线程并行计算耗时22秒。相同条件下SCF计算不用MPEC+COSX方法加速,计算耗时110秒, MPEC+COSX 大约加速了 5 倍。

对于含有重元素的体系的涉及梯度的计算(尤其是O1NumHess计算),默认的COSX格点可能不够大,导致结果有数值噪音,例如出现“假虚频”。此时可以增加COSX格点,并且/或者设严COSX格点的截断阈值。例如

$scf
...
cosxgrid
35 # first number: radial grid, second number: angular grid
   # default: 02
cosxthr # screening threshold, the higher, the tighter
6
$end