计算化学公社

标题: PySCF测评 [打印本页]
作者
Author: Accelerator    时间: 2018-8-1 05:38
标题: PySCF测评
本帖最后由 Accelerator 于 2018-8-4 23:00 编辑

    经大家指正并确认过后发现PySCF默认是全核并行的，因此实际上与Gaussian相比并无速度优势。

    PySCF是用Python写的量子化学程序。程序提供简单，轻量级的高效平台，用于量子化学计算和代码开发。当前版本可以在一个SMP节点上有效地操作中等尺寸体系（大约3000基函数，1000电子）。
（以上复制自量子化学网）

    PySCF的安装非常简单，运行pip install pyscf即可。用PySCF进行计算需要编写Python脚本。在这里以苯乙炔的单点能为例。
  1. <font face="宋体" size="2">import datetime
  2. startTime=datetime.datetime.now()
  3. from pyscf import gto
  4. mol=gto.Mole()
  5. mol.atom=' C                  3.23425500    0.00028800    0.00000500;\
  6. C                  2.02419900   -0.00014700    0.00000000;\
  7. H                  4.30046100   -0.00092300   -0.00003200;\
  8. C                  0.59427700   -0.00006400    0.00000300;\
  9. C                 -0.11989600    1.21305600   -0.00000200;\
  10. C                 -0.11999300   -1.21311200    0.00000200;\
  11. C                 -1.51250400    1.20860200    0.00000100;\
  12. H                  0.42865400    2.14989400   -0.00000300;\
  13. C                 -1.51261100   -1.20854100   -0.00000200;\
  14. H                  0.42847400   -2.15000300    0.00000200;\
  15. C                 -2.21304900    0.00005500   -0.00000100;\
  16. H                 -2.05292200    2.15133500    0.00000200;\
  17. H                 -2.05309200   -2.15123800   -0.00000200;\
  18. H                 -3.29963800    0.00010900   -0.00000100'
  19. mol.basis='def2tzvp'
  20. mol.max_memory=2000
  21. from pyscf import dft
  22. scf=dft.RKS(mol)
  23. scf.xc='pbe0'
  24. scf.kernel()
  25. endTime=datetime.datetime.now()
  26. 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/），找到了内置基组的列表。
(, 下载次数 Times of downloads: 166)
对cc和def2系列基组支持相当完备，cc支持到5-zeta，带弥散的只支持无止尽的八月。Def2系列包含一系列RI辅助基组；完整内置pc和pcseg系列。这些基组的输入形式比较自由，如cc-pvtz, ccpvtz等均可。
(, 下载次数 Times of downloads: 170)
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”模式会输出什么。
(, 下载次数 Times of downloads: 163)
每一次迭代都会输出这样一大串信息。由此可以看出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里是否有波函数稳定性分析的功能。
  1. from pyscf import gto, scf
  2. mol = gto.M(atom='O 0 0 0; O 0 0 1')
  3. mf = scf.UHF(mol)
  4. mf.verbose=4
  5. mf = scf.addons.dynamic_sz_(mf)
  6. mf.kernel()
  7. print('S^2 = %s, 2S+1 = %s' % mf.spin_square())
复制代码
PySCF的DFT有一个好玩的特性，即可以自定义泛函。对于范围分离泛函可以通过omega属性指定omega值。不过手册里没有提到双杂化泛函，由于其使用了Libxc的泛函库，去查了Libxc的手册，其中也没有双杂化相关内容。看起来是不能使用双杂化了。
  1. from pyscf import gto, dft
  2. mol = gto.M(atom='H  0  0  0; F  0.9  0  0', basis='6-31g')
  3. mf = dft.RKS(mol)
  4. mf.xc = 'HF*0.2 + .08*LDA + .72*B88, .81*LYP + .19*VWN'
  5. mf.kernel()
  6. mf.xc = 'HF*0.5 + .08*LDA + .42*B88, .81*LYP + .19*VWN'
  7. mf.kernel()
  8. mf.xc = 'HF*0.8 + .08*LDA + .12*B88, .81*LYP + .19*VWN'
  9. mf.kernel()
  10. mf.xc = 'HF'
  11. mf.kernel()
复制代码
PySCF的RI效率如何呢；在一开始的代码里SCF部分改成如下内容：
  1. <font face="宋体" size="2">from pymf import dft
  2. mf=dft.RKS(mol)
  3. mf.xc='pbe0'
  4. mf.verbose=4
  5. mf = mf.density_fit()
  6. 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)。写法如下：
  1. <font face="宋体" size="2">mol.basis='6-31G*'
  2. mol.max_memory=2000
  3. from pyscf import scf, mcscf
  4. mf = scf.RHF(mol).run()
  5. mc = mcscf.CASCI(mf, 4, 4)
  6. 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的讨论就转过来了。

作者
Author: sobereva    时间: 2018-8-1 06:01
Hessian误拼成了Heissian

关于"比Gaussian的默认收敛限还要严格一个数量级"，这其实得看具体怎么定义的，Gaussian里同时用密度矩阵和能量作为判断标准，默认标准下对密度矩阵的要求已经严过头了，达到收敛的时候，大多数情况下能量变化已经远小于1E-9了。

DFT积分格点的设置不同，可以导致耗时相差非常大。G09默认的fine比多数程序都要高。6-311+G(d,p)下两个程序的结果其实差异也不算很小，我没用过PySCF，但怀疑积分格点可能比G09默认的低不少。

另外，指令集的利用也明显影响耗时。G09没有支持AVX2的版本，PySCF的积分库应该充分利用了AVX2。

Gaussian当基组里有f角动量的时候，诸如def2-TZVP，耗时会比不带f的猛增很多，可能PySCF这点占了优势。

PySCF的一个特色应该算是可以方便地做DMRG。
作者
Author: zjxitcc    时间: 2018-8-1 09:47
我也感觉在高角动量时PySCF比高斯快很多。
PS1：在PySCF中使用自定义基组，令其基组数据完全与高斯一致，可使波函数方法（如HF, CASSCF）的电子能量与高斯在小数点后前6-8位一致（还依赖于收敛精度）。
PS2：PySCF有检测波函数稳定性的功能，例子如pyscf/examples/scf/17-stability.py。不过HF的收敛技巧和收敛性与高斯还无法相媲美。
PS3：自从用了PySCF的CASSCF之后，感觉飞上天，孙启明的黑科技二阶CASSCF啊
最后也推销一下自己参与的公众号“量子化学”，其中写了一系列安装PySCF和Block做DMRG及NEVPT2的帖子
作者
Author: ene    时间: 2018-8-1 17:03
本帖最后由 ene 于 2018-8-1 22:49 编辑

确实效率很高啊，看来值得使用
PS. 话说这么比，不怕banned by Gaussian吗
作者
Author: Warm_Cloud    时间: 2018-8-1 19:02
本帖最后由 Warm_Cloud 于 2018-8-1 19:07 编辑

很奇怪那个80s和318s 的，它是不是默认并行？
作者
Author: Accelerator    时间: 2018-8-1 21:49
zjxitcc 发表于 2018-8-1 09:47
我也感觉在高角动量时PySCF比高斯快很多。
PS1：在PySCF中使用自定义基组，令其基组数据完全与高斯一致， ...

不知道PySCF的CASSCF和molpro之类的专业软件比起来怎么样
作者
Author: Accelerator    时间: 2018-8-1 21:50
ene 发表于 2018-8-1 17:03
确实效率很高啊，看来值得使用
PS. 话说这么比，不怕baned by Gaussian吗

绝对会被ban的(x
被ban了就用ORCA 无所畏惧
作者
Author: Accelerator    时间: 2018-8-1 21:50
Warm_Cloud 发表于 2018-8-1 19:02
很奇怪那个80s和318s 的，它是不是默认并行？

这个确认过的，没有并行
作者
Author: ene    时间: 2018-8-1 22:36
Accelerator 发表于 2018-8-1 21:50
绝对会被ban的(x
被ban了就用ORCA 无所畏惧

大牛风范啊
作者
Author: wangxubo    时间: 2018-8-2 01:00
算casscf可以指定dmrg,fciqmc还有shci作为fcisolver,这样你的(54,54)就也能算了.
作者
Author: Warm_Cloud    时间: 2018-8-2 10:53
Accelerator 发表于 2018-8-1 21:50
这个确认过的，没有并行

我们自己测了一下，pyscf默认用并行，我们的电脑上耗时65s，高斯16耗时36s。你top一下再确认是不是有并行。
作者
Author: zjxitcc    时间: 2018-8-2 11:18
PySCF默认是并行的，如misc.num_threads(n=4)可以设定并行核数。当然其中有些步骤（如Einstein求和）是串行的，这时候num_threads不起作用。有时候有些任务时间太短了，看不到并行也是可能的。
作者
Author: 我本是个娃娃    时间: 2018-8-2 15:14
传一个Pyscf 1.4版本的说明书

(, 下载次数 Times of downloads: 454)
作者
Author: Accelerator    时间: 2018-8-4 22:59
感谢大家的更正。重新看了一下应该确实是默认全核并行的。这样的话PySCF实际上应该没有什么速度优势。
作者
Author: pyscf    时间: 2018-8-5 02:41
Accelerator 发表于 2018-8-4 22:59
感谢大家的更正。重新看了一下应该确实是默认全核并行的。这样的话PySCF实际上应该没有什么速度优势。

吹上了天 呵呵
本尊都没出来吹
作者
Author: k64_cc    时间: 2018-8-13 12:05
然而他们其实是拿这个研究postHF的……

PySCF的MP2效率倒是比Gaussian要高，不知道和ORCA比如何。
作者
Author: 函数与激情    时间: 2020-6-15 22:58
简单测试了下PySCF CASSCF的梯度和MOlcas的CASSCF梯度 发现差别有点大啊
作者
Author: zjxitcc    时间: 2020-6-15 23:53
函数与激情 发表于 2020-6-15 22:58
简单测试了下PySCF CASSCF的梯度和MOlcas的CASSCF梯度 发现差别有点大啊

先说说看能量差多少
作者
Author: 函数与激情    时间: 2020-6-16 12:29
zjxitcc 发表于 2020-6-15 23:53
先说说看能量差多少

由于我要计算一个相对较大分子的CASSCF梯度能量，传说Pyscf CASSCF效率不错，因此尝试一下下看看效果咋样 我先简单测试了个丙烯醛小分子之前优化出来的S1S0交叉结构处的能量和梯度，和MOLCAS和MOLPRO对比了一下。结果如图所示，三者的能量都满足S1和S0近简并的这个趋势，MOLCAS和MOLPRO的能量和梯度结果差别是很小的，pyscf S1处的梯度也还好吧，但怎么感觉和S0一样，不确定S0是不是我算错了呀，求大神指导一下。 (, 下载次数 Times of downloads: 107) (, 下载次数 Times of downloads: 79) (, 下载次数 Times of downloads: 92)
作者
Author: zjxitcc    时间: 2020-6-16 13:52
函数与激情 发表于 2020-6-16 12:29
由于我要计算一个相对较大分子的CASSCF梯度能量，传说Pyscf CASSCF效率不错，因此尝试一下下看看效果咋样 ...

没看懂你到底想算哪两个态。PySCF给出的两个根是最低的三重态和单重态，注意看自旋。在不同的程序中，有的时候可以限制所有根均为单重态做计算，有的时候也可以不限制。看实际问题。
而且就算是S0/S1，也最好用态平均。看不出你几个软件分别用的是特定态还是态平均。
作者
Author: 函数与激情    时间: 2020-6-16 14:40
本帖最后由 函数与激情 于 2020-6-16 14:41 编辑
zjxitcc 发表于 2020-6-16 13:52
没看懂你到底想算哪两个态。PySCF给出的两个根是最低的三重态和单重态，注意看自旋。在不同的程序中，有 ...

MOLCAS和MOLPRO的输入这里我没贴，两者用的都是S0,S1两态平均，PYSCF之前没接触过，在example里面找的输入，应该是没有使用态平均的，使用state_specific分别对于0和1。
作者
Author: Kirisame_Sakuya    时间: 2020-11-6 12:59
找到一份pyscf的中文文档，里面比较详细地介绍了pyscf的主要功能，并且利用pyscf的接口实现了XYG3解析梯度，对学习pyscf以及对自己亲手实现理论感兴趣的同学很有帮助。

https://py-xdh.readthedocs.io/zh_CN/latest/
作者
Author: jiangning198511    时间: 2020-11-6 13:55
Kirisame_Sakuya 发表于 2020-11-6 12:59
找到一份pyscf的中文文档，里面比较详细地介绍了pyscf的主要功能，并且利用pyscf的接口实现了XYG3解析梯度 ...

好东西 感谢
作者
Author: zkzk1234    时间: 2021-9-4 18:34
如果手动编译调用一些intel的数据库,不知道会不会快一点.
作者
Author: hebrewsnabla    时间: 2021-9-4 21:50
zkzk1234 发表于 2021-9-4 18:34
如果手动编译调用一些intel的数据库,不知道会不会快一点.

默认编译方式就是使用intel数学库的，如果cmake能找得到的话。

conda安装的话，也是链接intel数学库的。
作者
Author: j5888xm    时间: 2022-9-27 11:50
我本是个娃娃 发表于 2018-8-2 15:14
传一个Pyscf 1.4版本的说明书

大佬有2.0的新版pdf说明书吗，感谢
作者
Author: zjxitcc    时间: 2022-9-27 12:41
j5888xm 发表于 2022-9-27 11:50
大佬有2.0的新版pdf说明书吗，感谢

官方没有推出PDF，都是在线的https://pyscf.org/user.html#
作者
Author: j5888xm    时间: 2022-9-27 13:17
zjxitcc 发表于 2022-9-27 12:41
官方没有推出PDF，都是在线的https://pyscf.org/user.html#

好的，谢谢您



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