DIRAC：通过数值差分计算偶极矩

beefly

DIRAC可以在SCF/MCSCF/CI和闭壳层MP2级别计算标量、二分量、四分量相对论的偶极矩。对于其它理论方法，可以通过数值差分计算偶极矩。但是，由于DIRAC的总能量的定义不同于其它程序，计算步骤也稍有不同。

第一步：在输入文件中定义电场。例如，如果计算偶极矩z分量，需要在**HAMILTONIAN中插入

1. .OPERATOR
   
2. 'dipole'
   
3. DIAGONAL
   
4. 'ZDIPLEN'
   
5. COMFACTOR
   
6. +0.001

复制代码

或者简化形式

1. .OPERATOR
   
2. ZDIPLEN
   
3. COMFACTOR  +0.001

复制代码

第二步：提交计算，得到总能量E(+0.001)。

第三步：把第一步的电场强度改变符号，重新计算得到E(-0.001)。

第四步：数值差分，计算电子对偶极矩的贡献（μe）：μe = -ΔE/ΔF= -[E(+0.001)-E(-0.001)]/0.002

这里需要注意：对于使用了冻芯近似的post-HF计算， post-HF步骤的总能量不包含芯电子能量！此时必须分别计算HF对偶极矩的贡献（μHF）和关联能对偶极矩的贡献（μcorr）。

第五步：计算核排斥能对偶极矩的贡献（μnuc）。

如果在*.mol文件中没有指定对称性，程序会把分子转动、平移到标准方位。在输出文件最开始找到标准方位坐标，例如氟化氢分子的计算结果为

1. Centered and Rotated
   
2. --------------------
   
3.     1         0.00000000     0.00000000    -1.79453012       1
   
4.     9         0.00000000     0.00000000     0.09519602       1
   

复制代码

于是z分量的μnuc = 1*(-1.79453012) + 9*(0.09519602) = -0.93776594

如果在*.mol文件中指定了对称性，在输出文件最开始直接找到直角坐标，例如氟化氢分子的计算结果为

1.   Cartesian coordinates xyz format (angstrom)
   
2.   -------------------------------------------
   
3.     2
   
4. H_1    0.0000000000   0.0000000000   0.0000000000
   
5. F_1    0.0000000000   0.0000000000   1.0000000000
   

复制代码

于是z分量的μnuc = 1*(0.0000000000/c) + 9*(1.0000000000/c) = 17.00753521。其中c是Bohr到埃的换算因子0.529177，因为上面坐标的单位是埃。

注意：如果用了ECP，则要扣除芯电子个数。例如，F用了ECP（2个芯电子），则μnuc计算中的9*要改为7*。

第六步：计算总偶极矩：μ = μnuc + μe = μnuc + μHF [+ μcorr]，原子单位。换算为Debye需要再乘上2.541746。