XYG3泛函在常见软件中的使用方法(一)

2021-06-16 16:29:22 浏览数 (1)

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的结果如下

代码语言:javascript复制
=>"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.

0 人点赞