| 在加上电场时，DIRAC对总能量的定义不同于其它程序，因此数值差分计算步骤也稍有不同。 第一步：在输入文件中定义电场。如果计算偶极矩z分量，需要在**HAMILTONIAN中插入 .OPERATOR
'dipole'
DIAGONAL
'ZDIPLEN'
COMFACTOR
+0.001 |
或者简化形式 .OPERATOR
ZDIPLEN
COMFACTOR +0.001 |
第二步：提交计算，重新计算得到总能量E(+0.001)。 第三步：把第一步的电场强度改变符号，得到E(-0.001)。 注意对称性：(1) 若E(+0.001)与E(-0.001)对称等价，可以跳过此步，(2) 对于存在g/u对称性的体系，外电场会破坏g/u对称性。 第四步：数值差分，计算电子对偶极矩的贡献（μe）：μe = -ΔE/ΔF = -[E(+0.001)-E(-0.001)]/0.002（原子单位）。 特殊情况：使用冻芯近似的post-HF计算，post-HF总能量不包含芯电子能量！必须分别计算HF对偶极矩的贡献（μHF）和关联能对偶极矩的贡献（μcorr）。 第五步：计算核排斥能对偶极矩的贡献（μnuc）。分两种情况： ①如果在*.mol文件中没有指定对称性，程序会把分子转动、平移到标准方位。在输出文件最开始找到标准方位坐标，例如氟化氢分子的计算结果为 Centered and Rotated
--------------------
1 0.00000000 0.00000000 -1.79453012 1
9 0.00000000 0.00000000 0.09519602 1 |
于是z分量的μnuc = 1*(-1.79453012) + 9*(0.09519602) = -0.93776594 ②如果在*.mol文件中指定了对称性，在输出文件最开始直接找到直角坐标，例如氟化氢分子的计算结果为 Cartesian coordinates xyz format (angstrom)
-------------------------------------------
2
H_1 0.0000000000 0.0000000000 0.0000000000
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（最后一步对应冻芯post-HF的情况）。结果是原子单位，换算为Debye需要再乘上换算因子2.541746。 | 