Comparison between the ordering of Molecular Orbital coefficients between Gaussian and PySCF

Question

I am trying to compare the Molecular Orbital coefficients obtained from PySCF and G09 program. Below is the PySCF and Gaussian input and output. I am using 3-21G basis and RHF method for both the packages and the reference system in H3O+.

PySCF input

from pyscf import gto, scf, lib, cc, ao2mo, grad
from pyscf.gto import mole
mol = gto.M(
    verbose = 5,
    atom = [
           ['O' ,  0.074642 , -0.000000 , -0.026389],
           ['H' , -0.103145 ,  0.000000 ,  0.930145],
           ['H' ,  0.991917 ,  0.000000 , -0.350692],
           ['H' , -0.664854 , -0.000000 , -0.658618]],
    basis = '3-21g',charge=1,
cart=False
)
mf = scf.RHF(mol)
mf.conv_tol = 1e-13
mf.scf()
orbshape = mf.mo_coeff.shape
print('Orbital shape =',orbshape)
print("orbshape[1] is ", orbshape[1])
print(mf.mo_coeff)

PySCF output

[[ 9.83242599e-01 -2.36723812e-01 -3.60485158e-07  8.26285631e-08
   8.54201732e-18 -1.39400918e-01 -3.24924379e-07  4.27735532e-08
   7.40926127e-02 -5.10470621e-07 -2.15882620e-06  2.41476409e-18
  -5.78614062e-08  4.15554215e-07  7.81241453e-02]
 [ 9.72257975e-02  2.28321145e-01  4.37984919e-07 -3.38123813e-08
  -2.72406337e-18  7.77816722e-02  1.49262797e-07 -1.03367613e-07
  -9.97915817e-02  4.31272063e-07  3.64044620e-06 -1.81972141e-17
   9.28534624e-09 -1.68937209e-06 -1.68161970e+00]
 [-4.06623548e-02  7.36791646e-01  1.20392671e-06 -4.40830696e-07
  -5.31833825e-17  1.19408932e+00  3.50538923e-06  3.57091405e-07
  -6.69223322e-02  1.87873763e-06  1.99627276e-06  5.34167834e-17
   9.67059530e-07  2.32807183e-07  2.12095150e+00]
 [-1.36625076e-08 -3.18711971e-07  3.99875966e-01  1.75087384e-01
  -3.30919448e-16 -7.74299361e-07  1.65945660e-01 -2.49317047e-01
  -2.19084201e-06  7.50629546e-02 -1.99510623e-01  2.45189479e-16
  -3.12498671e-01 -1.00820556e+00  1.30919633e-06]
 [-2.20533621e-19  3.39181393e-18  3.32245076e-16  1.82208797e-16
   5.40820245e-01 -3.76638923e-17  2.52529300e-16  1.79190184e-16
  -3.58911407e-17 -1.86188339e-16  1.35022115e-16  1.01914696e+00
   6.37674753e-16  1.87552354e-16 -6.72338619e-18]
 [ 7.93352318e-09  3.52701876e-07 -1.75087393e-01  3.99876259e-01
  -1.46708126e-16 -4.71067753e-07  2.49317102e-01  1.65945561e-01
   2.10332685e-06  1.99503984e-01  7.50645771e-02  5.06756184e-16
  -1.00820614e+00  3.12498538e-01 -5.27814777e-07]
 [ 3.76724626e-08 -6.16668075e-07  4.07396435e-01  1.78380348e-01
  -2.65677104e-16 -1.38897953e-06  4.13014501e-01 -6.20514822e-01
  -9.57009076e-06  1.42830049e-01 -3.79622369e-01 -6.66840066e-16
   4.50878967e-01  1.45466938e+00 -3.93668672e-06]
 [ 9.64452306e-20  2.60830063e-19  4.13936439e-16  3.10077044e-16
   6.13588538e-01  6.65585084e-17 -1.99735894e-16 -4.35853887e-18
   7.78516109e-17  1.63777798e-16 -1.20469219e-16 -9.77065075e-01
  -5.86277720e-16 -1.23433162e-16  8.04573985e-18]
 [-2.39107820e-08  3.52553118e-07 -1.78380081e-01  4.07396534e-01
   5.72917657e-17 -1.06848036e-06  6.20513531e-01  4.13014793e-01
   4.28154714e-06  3.79622597e-01  1.42832608e-01 -6.40250772e-16
   1.45466023e+00 -4.50882758e-01  2.22748138e-06]
 [ 1.94248591e-03  9.23862276e-02 -1.28310475e-01  1.88977923e-01
   2.14775999e-16 -1.06720854e-01 -8.55193536e-02 -8.30944140e-02
   7.70539129e-01 -9.78680219e-01 -5.91504375e-01 -1.39238447e-16
  -2.47188467e-01  1.30059848e-01 -1.98379128e-01]
 [ 5.62456945e-03 -3.46785210e-04 -5.32958663e-02  7.84956157e-02
  -4.08253441e-17 -6.24400751e-01 -9.02544457e-01 -8.76971492e-01
  -4.28270887e-01  6.91196805e-01  4.17745391e-01  3.12184911e-16
  -5.00615506e-01  2.63385072e-01 -2.95165370e-01]
 [ 1.94245572e-03  9.23861587e-02  2.27814518e-01  1.66309493e-02
  -1.58478620e-16 -1.06722252e-01 -2.92034169e-02  1.15608088e-01
   7.70589446e-01 -2.29052135e-02  1.14327999e+00  2.02266614e-17
   1.09655908e-02 -2.79105624e-01 -1.98375611e-01]
 [ 5.62457575e-03 -3.46710029e-04  9.46264238e-02  6.90793578e-03
  -4.73459356e-17 -6.24403880e-01 -3.08207951e-01  1.22010974e+00
  -4.28309655e-01  1.61763449e-02 -8.07445431e-01  4.90902219e-16
   2.22085951e-02 -5.65236938e-01 -2.95165494e-01]
 [ 1.94248603e-03  9.23858781e-02 -9.95045334e-02 -2.05608780e-01
   1.65099861e-17 -1.06721571e-01  1.14722404e-01 -3.25123755e-02
   7.70550564e-01  1.00157189e+00 -5.51838435e-01  2.93965220e-16
   2.36222936e-01  1.49054040e-01 -1.98378874e-01]
 [ 5.62454904e-03 -3.46679780e-04 -4.13309319e-02 -8.54038412e-02
   1.17144303e-16 -6.24404361e-01  1.21074784e+00 -3.43140290e-01
  -4.28281297e-01 -7.07367252e-01  3.89732419e-01 -9.90418101e-16
   4.78406006e-01  3.01849336e-01 -2.95164069e-01]]

Gaussian input

# rhf/3-21g nosymm output=wfx pop=full

H3O+ run by HF/3-21G

1 1
O    0.074642   -0.000000   -0.026389
H   -0.103145    0.000000    0.930145
H    0.991917    0.000000   -0.350692
H   -0.664854   -0.000000   -0.658618

h3o+_hf.wfx

Gaussian output (only showing the Molecular orbital coefficients part)

     Molecular Orbital Coefficients:
                           1         2         3         4         5
                           O         O         O         O         O
     Eigenvalues --   -19.47020  -1.43467  -0.94637  -0.94637  -0.70968
   1 1   O  1S          0.98228  -0.23712   0.00000   0.00000   0.00000
   2        2S          0.10399   0.23148   0.00000   0.00000   0.00000
   3        2PX         0.00000   0.00000   0.40572   0.19088   0.00000
   4        2PY         0.00000   0.00000   0.00000   0.00000   0.54999
   5        2PZ         0.00000   0.00000  -0.19088   0.40572   0.00000
   6        3S         -0.04546   0.72117   0.00000   0.00000   0.00000
   7        3PX         0.00000   0.00000   0.37306   0.17551   0.00000
   8        3PY         0.00000   0.00000   0.00000   0.00000   0.60475
   9        3PZ         0.00000   0.00000  -0.17551   0.37306   0.00000
  10 2   H  1S          0.00202   0.09741  -0.13924   0.19360   0.00000
  11        2S          0.00620   0.00349  -0.06611   0.09192   0.00000
  12 3   H  1S          0.00202   0.09741   0.23729   0.02379   0.00000
  13        2S          0.00620   0.00349   0.11266   0.01129   0.00000
  14 4   H  1S          0.00202   0.09741  -0.09804  -0.21739   0.00000
  15        2S          0.00620   0.00349  -0.04655  -0.10321   0.00000
                           6         7         8         9        10
                           V         V         V         V         V
     Eigenvalues --    -0.24515  -0.13370  -0.13369   0.56303   0.63870
   1 1   O  1S         -0.14814   0.00000   0.00000   0.06752   0.00000
   2        2S          0.10957   0.00000   0.00000  -0.09812   0.00000
   3        2PX         0.00000  -0.09436  -0.35381   0.00000  -0.06009
   4        2PY         0.00000   0.00000   0.00000   0.00000   0.00000
   5        2PZ         0.00000  -0.35381   0.09436   0.00000  -0.22535
   6        3S          1.17146   0.00000  -0.00001   0.00409   0.00000
   7        3PX         0.00000  -0.18278  -0.68534   0.00000  -0.07354
   8        3PY         0.00000   0.00000   0.00000   0.00000   0.00000
   9        3PZ         0.00000  -0.68534   0.18278   0.00000  -0.27579
  10 2   H  1S         -0.14877   0.17199  -0.08189   0.76216   1.00935
  11        2S         -0.59347   1.05279  -0.50130  -0.46513  -0.83550
  12 3   H  1S         -0.14878  -0.01507   0.18990   0.76218  -0.08848
  13        2S         -0.59345  -0.09225   1.16239  -0.46515   0.07324
  14 4   H  1S         -0.14878  -0.15692  -0.10800   0.76216  -0.92087
  15        2S         -0.59346  -0.96054  -0.66108  -0.46513   0.76226
                          11        12        13        14        15
                           V         V         V         V         V
     Eigenvalues --     0.63871   1.00729   1.24182   1.24182   2.41635
   1 1   O  1S          0.00000   0.00000   0.00000   0.00000   0.07889
   2        2S          0.00000   0.00000   0.00000   0.00000  -1.67911
   3        2PX        -0.22535   0.00000  -0.24259  -0.99572   0.00000
   4        2PY         0.00000   1.01423   0.00000   0.00000   0.00000
   5        2PZ         0.06010   0.00000  -0.99572   0.24259   0.00000
   6        3S          0.00000   0.00000   0.00000   0.00000   2.13981
   7        3PX        -0.27578   0.00000   0.37297   1.53089   0.00000
   8        3PY         0.00000  -0.98256   0.00000   0.00000   0.00000
   9        3PZ         0.07355   0.00000   1.53088  -0.37298   0.00000
  10 2   H  1S         -0.48059   0.00000  -0.30132   0.13556  -0.20131
  11        2S          0.39781   0.00000  -0.52871   0.23785  -0.30309
  12 3   H  1S          1.11440   0.00000   0.03327  -0.32874  -0.20130
  13        2S         -0.92246   0.00000   0.05837  -0.57680  -0.30309
  14 4   H  1S         -0.63384   0.00000   0.26805   0.19319  -0.20131
  15        2S          0.52467   0.00000   0.47034   0.33895  -0.30309

I compared the output for both packages and found that, Gaussian prints the output as atom numbers followed by orbitals as in aufbau principle (probably, not tested for higher order atoms though). So, we obtain O (1S,2S,2P,3S,3P) H(1S,2S) and so on. Whereas PySCF builds them as atom numbers -> 1S,2S,3S etc. Then comes the P orbitals. Such that, O(1S,2S,3S,2P,3P) and then comes H(1S,2S). Right? Is there any generic rule for the different conventions (of different QM packages) to print MO coefficients?

Actually, I want to load the Molecular Orbitals from Gaussian calculation and reformat them in PySCF format. I want to do that to bypass costly PySCF calculations (I have asked this question on this previously).

Answer 1

Answer for basis function conventions:

What you've concluded is correct, but there are more rules for a general basis set (e.g. def2-TZVP, cc-pVTZ):

(1) The order of MO coefficients with angular momentum 5D/7F/9G/11H (and 6D/10F/15G/21H) of basis functions are usually different for different quantum chemistry packages, we need to permute these coefficients.

(2) For some MO coefficients, we need to multiply them by a normalized constant like sqrt(2), sqrt(3), etc.

(3) Although the basis set name 'cc-pVTZ' are specified in two packages, their primitive Gaussian exponents and contraction coefficients may be not the same. To solve this problem, we need to use exactly the same basis set data in two packages. The best way is to transfer MOs as well as basis set data from one package to another.

Conventions of basis functions in various quantum chemistry packages are recorded in the open source package MOKIT. For example, conventions of basis functions in PySCF/Gaussian and the related conversion relationship are coded in the file src/fch2py.f90; its reversed version is coded in the file src/py2fch.f90.

Also, you can find
fch2mkl/mkl2fch for Gaussian<->ORCA;
fch2inp/dat2fch for Gaussian<->GAMESS;
fch2com/xml2fch for Gaussian<->Molpro;
, etc. And py2xxx modules for PySCF->other packages. A lot of quantum chemistry packages are supported. It would be tedious to show all conversions here.

Answer for loading MOs of Gaussian in PySCF:

This can be easily achieved via using MOKIT. Assuming I have a Gaussian .fch(k) file h2o.fch right now, in which the RHF/cc-pVTZ calculation is performed. Here is a Python script:

from mokit.lib.gaussian import load_mol_from_fch, mo_fch2py
from mokit.lib.py2fch_direct import fchk
from pyscf import scf, mcscf

# load Cartesian coordinates, basis set data, charge, spin, etc
fchname = 'h2o.fch'        # converged RHF orbitals
fchname1 = 'h2o_cas44.fch' # converged CASSCF(4,4) orbitals
mol = load_mol_from_fch(fchname)

# generate the mf object, make all needed arrays allocated
mf = scf.RHF(mol)
mf.init_guess = '1e'
mf.max_cycle = 1
mf.verbose = 4
mf.kernel()

# load MO coefficients (work for RHF/ROHF/UHF-type wfn)
mf.mo_coeff = mo_fch2py(fchname)
mf.max_cycle = 10
dm = mf.make_rdm1() # generate density matrix
mf.kernel(dm0=dm)
# here the RHF is converged in 2 cycles with energy change less than 3e-8 a.u.

# if we have converged CASSCF orbitals, we can also perform CASSCF
mf.mo_coeff = mo_fch2py(fchname1)
mc = mcscf.CASSCF(mf, 4, 4)
mc.kernel()
# here the CASSCF is converged in 1 cycle

# We can also export PySCF MOs to Gaussian .fch file.
# fchk() will generate a .fch file from scratch.
mf.mo_coeff = mc.mo_coeff.copy()
fchk(mf, 'test.fch')

If you are working with 1-RDM, you can generate natural orbitals (NOs) using 1-RDM and then transfer NOs among various packages. By the way, the keyword nosymm int=nobasistransform are recommended to be written in the Gaussian input file for a better transferring of MOs.

Hope this will do some help.