计算化学公社

标题: DIRAC：通过数值差分计算偶极矩 [打印本页]
作者
Author: beefly    时间: 2016-12-27 11:53
标题: DIRAC：通过数值差分计算偶极矩
本帖最后由 beefly 于 2016-12-27 17:02 编辑

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。





欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3