[求助]pyscf与molpro结果比较问题

最爱喵星人

请问两个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?

zjxitcc

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

直接对比不同程序的分子轨道系数的话，要格外小心。不同程序的原子基重叠积分矩阵S未必一样；且默认用球谐型、Cartesian型哪种基函数未必一样。做个单参考CC，轨道当然还是正则轨道，你可以对比轨道图形是否一致。
当然，不同程序的STO-3G计算下得到的分子轨道系数一般是一样的这种可以直接对比。基组大了就复杂了。

wangxubo

楼上正解，pyscf里冻两个轨道就一样了

最爱喵星人

谢谢。确实是差了两个冻核轨道。

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

    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)

zjxitcc

最爱喵星人 发表于 2019-11-15 17:02
谢谢。确实是差了两个冻核轨道。

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

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