XYG3型泛函是一类重要的双杂化泛函,包括XYG3, XYGJ-OS, XYG7等。其中XYG3泛函由张颖、徐昕和W. A. Goddard III于2009年在PNAS杂志上发表。由于这类泛函使用了不同泛函来计算密度和能量(即能量泛函是非自洽的),使得用户往往不能简单地在常见程序中使用。目前,通过内置或非内置的形式可以做XYG3型泛函计算的程序包括但不限于:Q-Chem, Gaussian, PySCF, FHI-aims, ORCA, ...
下面我们介绍一下使用Gaussian和PySCF做XYG3型泛函计算的方法。
Gaussian
方法一:手动写Link1
下面的输入是冻核的XYG3计算
代码语言:javascript复制%chk=test.chk
#p B3LYP/6-311 G(3df,2p)
title
0 1
N 0.0 0.0 0.108827
H 0.0 0.947455 -0.253931
H 0.820520 -0.473727 -0.253931
H -0.820520 -0.473727 -0.253931
--link1--
%chk=test.chk
#p B2PLYP chkbasis Guess=Read Geom=AllCheck SCF(MaxCycle=-1)
IOp(3/76=1000008033,3/77=0210701967,3/78=0678906789,3/125=0321103211)
SCF(MaxCycle=-1)表示不做轨道优化,直接算能量。各个IOP的含义可在高斯官网http://gaussian.com/iops找到。注意这种写法只适用于G03, G09等。G16不支持这样的写法,因为G16要求MP2轨道必须是SCF收敛的,若不收敛会报SCF Error错误而强行终止程序。由于显示的原因可能行距略大,注意两行关键词之间无空行。正常结束后可在输出文件里找到结果
代码语言:javascript复制E(B2PLYP) = -0.56547860493069D 02
此即为XYG3电子能量。
方法二:使用xDH4Gau
张颖等最近开源的xDH4Gau程序支持更多的xDH型泛函(如XYGJ-OS, XYG7等)的单点计算,可调用G03、G09、G16,也支持使用高斯的PCM等功能。程序地址:
代码语言:javascript复制https://github.com/igor-1982/xDH4Gau
安装、使用方法见该项目的README。例如,去掉Tests/Test001.gjf
中的full之后,做一个冻核的XYG3的结果如下
=>"XYG3" is choosen for the question
==----------------------------------------------------------------------------==
E(B3LYP) = -56.58684505 A.U. E(XYG3) = -56.54786054 A.U.
==----------------------------------------------------------------------------==
=>XYG3 belongs to the family of xDH@B3LYP
=>E(XYGJ-OS) = -56.42918257 A.U.
=>E(revXYG3) = -56.55109989 A.U.
=>E(XYG5) = -56.45854916 A.U.
=>E(XYG6) = -56.44674223 A.U.
=>E(XYG7) = -56.25195145 A.U.
=>Job Type :: Single-Point Calculation
不论用户选择哪种泛函,程序都会输出所有可用的XYG3型泛函的结果。
使用PySCF
尽管PySCF没有内置任何的双杂化泛函,但是只要熟悉双杂化泛函的逻辑,就能利用PySCF写出几行代码的XYG3运行脚本,这在PySCF的一个issue中有详尽的讨论:
代码语言:javascript复制https://github.com/pyscf/pyscf/issues/911
例如一个非冻核的XYG3计算:
代码语言:javascript复制from pyscf import gto, dft, mp
def get_energy_xDH(mol: gto.Mole, xc_scf: str, c_pt2: float, xc_dh: str or None=None):
# SCF part calculation
mf_scf = dft.KS(mol, xc=xc_scf)
mf_scf.conv_tol_grad = 1e-10
mf_scf.grids.atom_grid = (99, 590)
mf_scf.run()
if not mf_scf.converged:
raise ValueError("SCF not converged.")
e_nscf = mf_scf.e_tot
# Non-consistent part calculation
if xc_dh is not None:
mf_nscf = dft.KS(mol, xc=xc_dh)
mf_nscf.grids.atom_grid = (99, 590)
e_nscf = mf_nscf.energy_tot(dm=mf_scf.make_rdm1())
# PT2 contribution
mf_pt2 = mp.MP2(mf_scf).run()
e_pt2 = c_pt2 * mf_pt2.e_corr
return e_nscf e_pt2
mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build()
# XYG3 energy : -76.2910536679
e_xyg3 = get_energy_xDH(mol, "B3LYPg", 0.3211, "0.8033*HF - 0.0140*LDA 0.2107*B88, 0.6789*LYP")
print(e_xyg3)
我们也可以利用现成的Py_xDH程序来完成这样的计算,项目地址:
代码语言:javascript复制https://github.com/ajz34/Py_xDH
Py_xDH是目前唯一的可做XYG3梯度、Hessian和极化率的公开程序。该程序的文档
代码语言:javascript复制https://py-xdh.readthedocs.io/zh_CN/latest/qcbasic/basic_xyg3.html
提供了对XYG3定义的详细描述。由于该程序对于开发者以外的用户来说可能难以使用,从xDH4Gau到Py_xDH的接口程序正在开发中。
相关文献 [1] Zhang, Y.; Xu, X.; W. A. Goddard. Doubly hybrid density functional for accurate descriptions of nonbond interactions, thermochemistry, and thermochemical kinetics. PNAS 2009, 106, 4963. [2] 张颖, 徐昕. 新一代密度泛函方法XYG3[J]. 化学进展, 2012, 24(006):1023-1037. [3] Zhang, I. Y.; Xu, X. Exploring the Limits of the XYG3-Type Doubly Hybrid Approximations for the Main-Group Chemistry: The xDH@B3LYP Model. JPCL 2021, 12, 2638.