本帖最后由 Accelerator 于 2018-8-4 23:00 编辑
经大家指正并确认过后发现PySCF默认是全核并行的,因此实际上与Gaussian相比并无速度优势。
PySCF是用Python写的量子化学程序。程序提供简单,轻量级的高效平台,用于量子化学计算和代码开发。当前版本可以在一个SMP节点上有效地操作中等尺寸体系(大约3000基函数,1000电子)。 (以上复制自量子化学网)
PySCF的安装非常简单,运行pip install pyscf即可。用PySCF进行计算需要编写Python脚本。在这里以苯乙炔的单点能为例。- <font face="宋体" size="2">import datetime
- startTime=datetime.datetime.now()
- from pyscf import gto
- mol=gto.Mole()
- mol.atom=' C 3.23425500 0.00028800 0.00000500;\
- C 2.02419900 -0.00014700 0.00000000;\
- H 4.30046100 -0.00092300 -0.00003200;\
- C 0.59427700 -0.00006400 0.00000300;\
- C -0.11989600 1.21305600 -0.00000200;\
- C -0.11999300 -1.21311200 0.00000200;\
- C -1.51250400 1.20860200 0.00000100;\
- H 0.42865400 2.14989400 -0.00000300;\
- C -1.51261100 -1.20854100 -0.00000200;\
- H 0.42847400 -2.15000300 0.00000200;\
- C -2.21304900 0.00005500 -0.00000100;\
- H -2.05292200 2.15133500 0.00000200;\
- H -2.05309200 -2.15123800 -0.00000200;\
- H -3.29963800 0.00010900 -0.00000100'
- mol.basis='def2tzvp'
- mol.max_memory=2000
- from pyscf import dft
- scf=dft.RKS(mol)
- scf.xc='pbe0'
- scf.kernel()
- endTime=datetime.datetime.now()
- print (endTime-startTime).seconds</font>
复制代码 首先通过gto定义分子构型。gto的属性包含了基组(basis),电荷(charge),自旋多重度(spin),同位素(mass)等,定义方式非常自然。例如基组既可以用上述字符串定义,如果希望使用混合基组,也可以用数组定义。以下是官网给出的例子。
mol.symmetry = 1mol.charge = 1mol.spin = 1mol.nucmod = {'O1': 1}mol.mass = {'O1': 18, 'H': 2} PySCF的官方文档上写分子构型中不同原子可以通过分号或换行分隔,所以如果直接读取xyz文件的话应该可以不作处理。不过由于字符串中间直接换行是不合语法的,像现在这样手动输入坐标的话就还是要在换行处补加斜杠了,因此也顺便加上了分号。PySCF对笛卡尔坐标和内坐标均支持。
PySCF官网并没有列出其内置的基组。Rita一开始尝试在basis中指定6-311+G(d,p)报错后在报错信息中找到了内置基组的存放位置(/usr/local/lib/python2.7/dist-packages/pyscf/gto/basis/),找到了内置基组的列表。
对cc和def2系列基组支持相当完备,cc支持到5-zeta,带弥散的只支持无止尽的八月。Def2系列包含一系列RI辅助基组;完整内置pc和pcseg系列。这些基组的输入形式比较自由,如cc-pvtz, ccpvtz等均可。
Pople系列基组的文件名称就显得有些诡异。经过尝试,输入6-31G(d)可以识别,但6-311+G(d,p)就不行,要写作6-311+G**. 执行上述文件输出到test1.log里,可以看到计算得到的能量,用时80 s. 默认的输出详细度下只输出了能量的数值。 作为能量和效率的对比,用Gaussian进行相同的计算(并设置使用的内存同为2000 MB)。G09给出的能量是-308.124621981,与PySCF基本相同,用时为318 s。 为了一探Gaussian为什么比PySCF慢这么多,把PySCF的输出详细程度调到最大(scf.verbose=9,一般设到4或5即可得到需要的信息),同时也看一看PySCF的“noisy”模式会输出什么。 每一次迭代都会输出这样一大串信息。由此可以看出SCF收敛限为1E-9,比Gaussian的默认收敛限还要严格一个数量级;整个任务经过11次迭代,比Gaussian少4次,不过也不足以解释巨大的时间差。因此PySCF的速度优势应当是在算法的某处。所以即使Gaussian的双电子积分几乎是令人绝望的优秀,但在其他地方的效率还是有很多可以优化之处的。顺带一提PySCF的DFT默认没有开RI(可以通过with_df.auxbasis属性设置),如果再开RI应该还会更快。 尝试了def2基组之后,再来看一看PySCF对Gaussian最擅长的Pople基组的表现。6-311+G(d,p)在PySCF里要写作6-311+G**。本次计算用时57 s,能量-308.090572698915;Gaussian用时43 s,能量-308.090450328。Gaussian对Pople基组的效率还是相当好的。 接下来测试UDFT。首先是苯乙炔的三重态,用PBE0/def2-SVP算单点。注意这里的mol.spin指定的是alpha电子减去beta电子,因此等于自旋多重度-1.这里输出详细度指定为4,就只输出了每一轮迭代的能量。计算在22 s完成,能量为-307.661565969255。Gaussian在相同水平下的计算用时51 s,能量-307.665874496。在官方手册上波函数稳定性分析部分神奇地留空了。在下方给出了一个自动调节自旋态的功能,在这个例子中虽然没有指定氧气的自旋态,但如果使用UHF,通过addons的设置可以自动收敛到三重态基态。很想试验一下这个功能,但发现似乎只能用于UHF,不支持将UKS的mf对象作为参数(之前的代码中把scf用作了对象的名称,但现在如果import scf的话就只能用其他名称了)。由于手册里相应部分留空,很好奇PySCF里是否有波函数稳定性分析的功能。 - from pyscf import gto, scf
- mol = gto.M(atom='O 0 0 0; O 0 0 1')
- mf = scf.UHF(mol)
- mf.verbose=4
- mf = scf.addons.dynamic_sz_(mf)
- mf.kernel()
- print('S^2 = %s, 2S+1 = %s' % mf.spin_square())
复制代码PySCF的DFT有一个好玩的特性,即可以自定义泛函。对于范围分离泛函可以通过omega属性指定omega值。不过手册里没有提到双杂化泛函,由于其使用了Libxc的泛函库,去查了Libxc的手册,其中也没有双杂化相关内容。看起来是不能使用双杂化了。 - from pyscf import gto, dft
- mol = gto.M(atom='H 0 0 0; F 0.9 0 0', basis='6-31g')
- mf = dft.RKS(mol)
- mf.xc = 'HF*0.2 + .08*LDA + .72*B88, .81*LYP + .19*VWN'
- mf.kernel()
- mf.xc = 'HF*0.5 + .08*LDA + .42*B88, .81*LYP + .19*VWN'
- mf.kernel()
- mf.xc = 'HF*0.8 + .08*LDA + .12*B88, .81*LYP + .19*VWN'
- mf.kernel()
- mf.xc = 'HF'
- mf.kernel()
复制代码PySCF的RI效率如何呢;在一开始的代码里SCF部分改成如下内容: - <font face="宋体" size="2">from pymf import dft
- mf=dft.RKS(mol)
- mf.xc='pbe0'
- mf.verbose=4
- mf = mf.density_fit()
- mf.kernel()</font>
复制代码以默认的辅助基组进行RI.用时只降低了4 s。在mf =mf.density_fit()后面写mf.with_df.auxbasis = 'def2-tzvp-ri'指定辅助基组,耗时也只从80 s降低到了74 s。不过不加RI本身就已经这么快了,这样的成绩也算是可以接受。 然后是PySCF的CASSCF功能。仍然以苯乙炔为例,考虑到CASSCF一般比较慢,所以基组小至6-31G(d),活性空间取(4,4)。写法如下: - <font face="宋体" size="2">mol.basis='6-31G*'
- mol.max_memory=2000
- from pyscf import scf, mcscf
- mf = scf.RHF(mol).run()
- mc = mcscf.CASCI(mf, 4, 4)
- mc.kernel()</font>
复制代码令人震惊的是仅仅3 s计算就结束了。得到能量为-306.401162634104。索性将活性空间扩大到(10,10),计算用时4 s, 能量为-306.42444538471。接下来将活性空间扩大到(54,54)(PySCF也有FCI功能),10分钟没有算完,放弃。此外PySCF也可以手动指定纳入活性空间的轨道标号。 相较之下Gaussian在(4,4)的活性空间下用时119 s,收敛失败。xswl 除此之外PySCF里还有一些有用的特性: 可以自定义哈密尔顿 可以生成任意原子轨道积分;此外PySCF还提供了libcint积分库的接口可以自己退(划掉)自己用 可以做轨道定域化 可以算周期性体系(手册上写了很多,不知道这个功能怎么样) 可以提取解析梯度(手册上解析梯度只介绍了HF,而Heissian虽然有这一条目但也是空的。尝试写grad.RKS但报错了)
PySCF不支持的功能 隐式溶剂化 既然梯度和Hessian那么有限应该不能用来做几何结构优化和频率计算 之前提到过的双杂化泛函 耦合簇只支持CCSD和CCSD(T),不支持DLPNO(但支持F12)
由此来看,Pyscf如果要与主流量化软件抗衡还远远不够。不过由于其开源以及Python易于扩展的特性,很高的计算效率,以及可以调泛函调积分等功能,作为新的泛函等的开发工具,或是为一些需求不超出其功能范围的程序提供快速计算,无疑是十分有吸引力的。
本文首先发表在了自己的公众号“Rita是美少女”上;鉴于论坛上似乎还没有有关PySCF的讨论就转过来了。
|