Multireference Tutorial

With pyAutoMR

Current strategy of pyAutoMR:

If RHF is stable: use RHF -> PM LMO -> CASSCF;

Otherwise, use UHF -> UNO -> CASSCF.

“RHF stable” means E(RHF) = E(UHF).

UHF case

from pyscf import lib
import guess, autocas, cidump

lib.num_threads(4)

xyz = 'N 0.0 0.0 0.0; N  0.0 0.0 1.9' #sys.argv[1]
bas = 'cc-pvdz'

mf = guess.from_frag(xyz, bas, [[0],[1]], [0,0], [3,-3], cycle=50)
guess.check_stab(mf)

mf2 = autocas.cas(mf)
cidump.dump(mf2)

For a stretched N2, we need to do UHF first.

Note

Simple UHF guess will fall into RHF, and even guess=mix maybe doesn’t work. A practical way is using fragment guess, i.e. guess=fragment in Gaussian, or guess.from_frag in pyAutoMR.

Here, [[0],[1]] means the 0th atom belongs to frag0, and 1st atom belongs to frag1. [0,0] means charges of the 2 frags are 0,0 respectively. [3,-3] means the spin of the 2 frags. Note that spin equals 2S, or n_a - n_b in PySCF, not 2S+1. Also, we need to ensure that the total number of alpha/beta electrons equals that in total molecule, so the spin parameter cannot be [3,3]. You can see output like

**** generating fragment guess ****
fragments: [('N', [0.0, 0.0, 0.0])] [('N', [0.0, 0.0, 1.9])]
converged SCF energy = -54.3911145621998  <S^2> = 3.7540306  2S+1 = 4.0020148
converged SCF energy = -54.3911145621998  <S^2> = 3.7540306  2S+1 = 4.0020148
       na   nb
atom1   5    2
atom2   2    5

Then, check_stab will do something like stable=opt in Gaussian, and autocas.cas will do UHF MO -> UNO, and estimate the size of active space. Output file contains a clear illustration of the UNOs

UNO ON: [ 2.        2.        1.999103  1.997838  1.550306  1.174501  1.174501  0.825499  0.825499  0.449694  0.002162  0.000897  0.        0.        0.
  0.        0.        0.        0.        0.       -0.       -0.       -0.       -0.       -0.       -0.       -0.       -0.      ]
UNO in active space
             #0      #1      #2      #3      #4      #5     
  0 N 3s                                               0.184
  0 N 2px              0.154  -0.432  -0.487  -0.139        
  0 N 2py              0.432   0.154  -0.139   0.487        
  0 N 2pz     -0.442                                   0.565
  0 N 3px              0.106  -0.295  -0.323                
  0 N 3py              0.295   0.106           0.323        
  0 N 3pz     -0.295                                   0.390
  1 N 3s                                              -0.184
  1 N 2px              0.154  -0.432   0.487   0.139        
  1 N 2py              0.432   0.154   0.139  -0.487        
  1 N 2pz      0.442                                   0.565
  1 N 3px              0.106  -0.295   0.323                
  1 N 3py              0.295   0.106          -0.323        
  1 N 3pz      0.295                                   0.390

pyAutoMR will pick the orbitals having occupation number (ON) bewteen 1.98-0.02. And it’s clear that the six orbitals are: simga bonding, pi bonding, pi bonding, pi anti-bonding, pi anti-bonding, and sigma anti-bonding.

Without pyAutoMR, we need to manually view and pick these orbitals, but now it’s automatic.