计算化学公社

标题: [求助]pyscf与molpro结果比较问题 [打印本页]
作者
Author: 最爱喵星人    时间: 2019-11-15 11:20
标题: [求助]pyscf与molpro结果比较问题
本帖最后由 最爱喵星人 于 2019-11-15 11:21 编辑

请问两个pyscf与molpro结果比较的问题：


pyscf输入文件：


import numpy
from pyscf import gto, scf,symm,cc,tools
from pyscf.cc import ccsd_t
mol=gto.Mole()
#unit in angstrom
mol.atom='''N 0 0 0
          N 0 0 1'''

mol.basis='def2tzvp'
mol.max_memory=4000
mol.charge=0
mol.spin=0  #molpro convention
mol.output='/home/test/pyscf/N2.out'
mol.verbose=4
mol.symmetry=0
mol.build()
hf=scf.RHF(mol).run()
#ccsd
ccsd=cc.CCSD(hf).run()
ccsdt=ccsd.ccsd_t()

#print orbital
hf.analyze(hf)
--------------------------------------------------
molpro 输入文件：
***,test
print,basis,orbitals   ! this is optional: print the basis set and
                       ! the occupied orbitals
geometry={             ! define the nuclear coordinates
N 0 0 0
N 0 0 1
}

basis=def2-tzvp              ! Select basis set
hf                     ! Invoke Hartree-Fock program
ccsd
---


pyscf HF 能量 ：-108.973429768344
molpro HF 能量： -108.97342977


但是
pcscf CCSD(T) 能量：
CCSD converged
E(CCSD) = -109.3649761375377  E_corr = -0.3915463691937363
CCSD(T) correction = -0.0152932556653192


molpro CCSD(T) 能量：

  CCSD total energy                  -109.326600584694
  CCSD[T] energy                     -109.342087953176
  CCSD-T energy                      -109.341347398788
!CCSD(T) total energy               -109.341596769085



请问我的pcscf哪里输错了导致二者ccsd和ccsd(t)能量不一样？


另外在打印出的分子轨道中：
pyscf：
               #1        #2        #3        #4        #5      
  0 N 1s       0.70709   0.70710   0.00218   0.00227   0.00000
  0 N 2s       0.00417   0.00118  -0.57335  -0.50802  -0.00000
  0 N 3s       0.00216   0.00022   0.06529  -0.00713   0.00000
  0 N 4s      -0.00030  -0.00064  -0.00418   0.00656  -0.00000
  0 N 5s      -0.00051  -0.00068   0.01057  -0.00133   0.00000
  0 N 2px      0.00000  -0.00000  -0.00000   0.00000   0.51493
  0 N 2py      0.00000  -0.00000   0.00000  -0.00000  -0.47977
  0 N 2pz      0.00145   0.00202   0.39832  -0.49163   0.00000
  0 N 3px     -0.00000   0.00000   0.00000  -0.00000  -0.01985
  0 N 3py     -0.00000   0.00000  -0.00000   0.00000   0.01849
  0 N 3pz     -0.00063   0.00028  -0.05591  -0.00493  -0.00000
  0 N 4px      0.00000  -0.00000   0.00000   0.00000  -0.00639
  0 N 4py      0.00000  -0.00000   0.00000  -0.00000   0.00596
  0 N 4pz      0.00039   0.00038   0.01201  -0.00445   0.00000
  0 N 3dxy    -0.00000  -0.00000   0.00000   0.00000  -0.00000
  0 N 3dyz    -0.00000  -0.00000  -0.00000   0.00000  -0.04100
  0 N 3dz^2    0.00120   0.00082   0.06670  -0.00740  -0.00000
  0 N 3dxz     0.00000  -0.00000  -0.00000   0.00000   0.04401
0 N 3dx2-y2   0.00000   0.00000   0.00000   0.00000  -0.00000
  0 N 4dxy    -0.00000  -0.00000  -0.00000  -0.00000   0.00000
  0 N 4dyz     0.00000  -0.00000  -0.00000   0.00000   0.00462
  0 N 4dz^2   -0.00024  -0.00013   0.00318  -0.00091  -0.00000
  0 N 4dxz     0.00000  -0.00000  -0.00000   0.00000  -0.00496
0 N 4dx2-y2  -0.00000  -0.00000  -0.00000  -0.00000   0.00000
  0 N 4fy^3    0.00000   0.00000   0.00000   0.00000  -0.00000
  0 N 4fxyz    0.00000   0.00000   0.00000  -0.00000   0.00000
  0 N 4fyz^2  -0.00000  -0.00000  -0.00000  -0.00000  -0.00952
  0 N 4fz^3    0.00082  -0.00010   0.02200   0.00284  -0.00000
  0 N 4fxz^2   0.00000  -0.00000  -0.00000   0.00000   0.01022
  0 N 4fzx^2   0.00000   0.00000  -0.00000   0.00000  -0.00000
  0 N 4fx^3   -0.00000  -0.00000  -0.00000   0.00000  -0.00000
  1 N 1s       0.70709  -0.70710   0.00218  -0.00227  -0.00000
  1 N 2s       0.00417  -0.00118  -0.57335   0.50802  -0.00000
  1 N 3s       0.00216  -0.00022   0.06529   0.00713  -0.00000
  1 N 4s      -0.00030   0.00064  -0.00418  -0.00656   0.00000
  1 N 5s      -0.00051   0.00068   0.01057   0.00133  -0.00000
  1 N 2px     -0.00000   0.00000  -0.00000   0.00000   0.51493
  1 N 2py      0.00000  -0.00000   0.00000   0.00000  -0.47977
  1 N 2pz     -0.00145   0.00202  -0.39832  -0.49163   0.00000
  1 N 3px      0.00000   0.00000   0.00000  -0.00000  -0.01985
  1 N 3py     -0.00000   0.00000  -0.00000   0.00000   0.01849
  1 N 3pz      0.00063   0.00028   0.05591  -0.00493  -0.00000
  1 N 4px     -0.00000   0.00000   0.00000  -0.00000  -0.00639
  1 N 4py      0.00000  -0.00000  -0.00000   0.00000   0.00596
  1 N 4pz     -0.00039   0.00038  -0.01201  -0.00445   0.00000
  1 N 3dxy     0.00000  -0.00000   0.00000  -0.00000   0.00000
  1 N 3dyz     0.00000  -0.00000  -0.00000   0.00000   0.04100

---------------------------------

molpro

                                                
   1 1s
  1 1s
  1 1s
  1 1s
  1 1s
  1 2pz
  1 2pz
  1 2pz
1 3d0
  1 3d2+
    1 3d0
   1 3d2+
   1 4f0
   1 4f2+
   2 1s
   2 1s
   2 1s
   2 1s
  2 1s
   2 2pz
    2 2pz
   2 2pz
   2 3d0
   2 3d2+
   2 3d0
   2 3d2+
   2 4f0
   2 4f2+












0.31718
0.44529
0.01508
-0.00441
0.00111
0.00214
-0.00086
0.00022
0.00076
0.00000
-0.00064
0.00000
0.00010
0.00000
0.31718
0.44529
0.01508
-0.00441
0.00111
-0.00214
0.00086
-0.00022
0.00076
0.00000
-0.00064
0.00000
-0.00010
0.00000












0.31753
0.44583
0.01437
0.00413
0.00759
0.00175
0.00586
0.00234
0.00081
0.00000
0.00230
0.00000
0.00046
0.00000
-0.31753
-0.44583
-0.01437
-0.00413
-0.00759
0.00175
0.00586
0.00234
-0.00081
0.00000
-0.00230
0.00000
0.00046
0.00000












-0.07801
-0.18003
0.26714
0.22675
0.03281
0.19228
0.09170
0.00692
0.02178
0.00000
0.01835
0.00000
0.00455
0.00000
-0.07801
-0.18003
0.26714
0.22675
0.03281
-0.19228
-0.09170
-0.00692
0.02178
0.00000
0.01835
0.00000
-0.00455
0.00000












0.06641
0.15720
-0.24259
-0.35388
-0.19895
0.16628
0.14779
0.05285
0.00438
0.00000
0.00880
0.00000
0.00064
0.00000
-0.06641
-0.15720
0.24259
0.35388
0.19895
0.16628
0.14779
0.05285
-0.00438
0.00000
-0.00880
0.00000
0.00064
0.00000












0.02299
0.05519
-0.07323
-0.19927
-0.18754
0.34327
0.23688
0.09709
0.01882
0.00000
0.01039
0.00000
0.00247
0.00000
0.02299
0.05519
-0.07323
-0.19927
-0.18754
-0.34327
-0.23688
-0.09709
0.01882
0.00000
0.01039
0.00000
-0.00247
0.00000







可以看到打印出分子轨道信息不一样。pyscf得到的分子轨道似乎进行局域化了。是这样么？如果是，怎么得到canonical orbital?




作者
Author: zjxitcc    时间: 2019-11-15 13:22
本帖最后由 zjxitcc 于 2019-11-15 13:24 编辑

做MP2, CC计算，写输入文件首先要想到的就是是否冻核（芯）、要冻结几个芯轨道。这个能量差异来自于有无冻芯。

直接对比不同程序的分子轨道系数的话，要格外小心。不同程序的原子基重叠积分矩阵S未必一样；且默认用球谐型、Cartesian型哪种基函数未必一样。做个单参考CC，轨道当然还是正则轨道，你可以对比轨道图形是否一致。
当然，不同程序的STO-3G计算下得到的分子轨道系数一般是一样的这种可以直接对比。基组大了就复杂了。
作者
Author: wangxubo    时间: 2019-11-15 13:49
楼上正解，pyscf里冻两个轨道就一样了
作者
Author: 最爱喵星人    时间: 2019-11-15 17:02
谢谢。确实是差了两个冻核轨道。

关于第二个问题，看起来是做了一个变换。

    if verbose >= logger.DEBUG:
        label = mf.mol.ao_labels()
        if with_meta_lowdin:
            log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) **')
            orth_coeff = orth.orth_ao(mf.mol, 'meta_lowdin', s=ovlp_ao)
            c = reduce(numpy.dot, (orth_coeff.conj().T, ovlp_ao, mo_coeff))
        else:
            log.debug(' ** MO coefficients (expansion on AOs) **')
            c = mo_coeff
        dump_mat.dump_rec(mf.stdout, c, label, start=MO_BASE, **kwargs)
    dm = mf.make_rdm1(mo_coeff, mo_occ)

作者
Author: zjxitcc    时间: 2019-11-15 18:01
最爱喵星人 发表于 2019-11-15 17:02
谢谢。确实是差了两个冻核轨道。

关于第二个问题，看起来是做了一个变换。

可以试试用tools.dump_mat.dump_mo，而非analyze



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