ROHF with contrained configuration

[Alexandre Voute]

Hello everyone,

[[attached output files are referenced in square brackets below]]

I am using Molpro 2020.2

I am working with the system [N2H + Rb](+), which is a singly-charged, doublet system. Here I keep the two moieties about 30Å far from each other. I use a large-core ECP (+CPP) for Rb so as to have to deal with only one electron in its 5s orbital.

I know from previous work that the lowest-energy configuration is the one where the unpaired electron is on the N2H moiety, i.e. we have N2H + Rb(+).

1) However, when I do a simple ROHF...

{rhf; wf,15,1,1}

...the optimization places that electron on the Rb, although the LUMO of the N2H(+) moiety (orbital #8.1) is lower in energy that the 5s orbital (#7.1) of Rb. [[n2hrb+_rohf_test.out]]

 ELECTRON ORBITALS
 =================

  Orbital  Occupation    Energy     Coefficients
...

-0.28581     1 2pz   0.75801     2 2pz  -0.39162
    6.1     2.00000    -0.97428     1 2px   0.59155     2 2px   0.71080
    7.1     1.00000    -0.17163     4 1s    0.63053     4 1s    0.39469
    8.1     0.00000    -0.19234     1 2px   0.69906     1 2px   0.30779     2 2px  -0.56673     2 2px  -0.30991
   ...

I thus try to play with the OPEN and CLOSED cards so as to enforce the occupation of that N2H(+) LUMO by the unpaired electron.

(I work in the Cs point group, so there are only 2 irreps)

2) With {rhf; wf, 15,1,1; open, 8.1 }, it gets worse, as the bonding pi_y orbital (#1.2) of N2H+ is emptied and the unpaired electron is still on the Rb [[n2hrb+_rohf_open_test.out]]

 ELECTRON ORBITALS
 =================

  Orbital  Occupation    Energy     Coefficients
    1.1     2.00000   -16.12125     2 1s    0.99866
    2.1     2.00000   -15.97859     1 1s    0.99939
    3.1     2.00000    -1.87753     1 1s    0.57448     1 2pz   0.35579     2 1s    0.69865     2 2pz  -0.35685
    4.1     2.00000    -1.34219     1 1s   -0.37425     2 1s    0.55604     2 2pz   0.61910     3 1s    0.61944     3 1s   -0.27774
    5.1     2.00000    -1.08818     1 2px   0.58032     2 2px   0.80783
    6.1     2.00000    -0.97716     1 1s   -0.67905     1 1s   -0.37478     1 2pz   0.67808     2 2pz  -0.38941
    7.1     2.00000    -0.62538     1 2px   0.86147     2 2px  -0.68789
    8.1     1.00000    -0.17174     4 1s    0.63053     4 1s    0.39469
    9.1     0.00000    -0.15307     1 1s   -0.69336     2 1s    0.29958     2 1s    0.41141     2 1s    0.81955     2 2pz  -0.37793
                                    3 1s   -0.25507     3 1s   -0.31823     3 1s   -0.66899
...
    1.2     0.00000    -0.54343     1 2py   0.46619     2 2py   0.65703
...

3) Thus i need to specify the doubly-occupied orbitals. However, with {rhf; wf, 15,1,1; close, 6, 1; open, 8.1 } an extra electron appeears in the occupation! [[n2hrb+_rohf_open-close_test.out]]

 PROGRAM * RHF-SCF (OPEN SHELL)       Authors: W. Meyer, H.-J. Werner


 NUMBER OF ELECTRONS:       8+    7-    SPACE SYMMETRY=1    SPIN SYMMETRY: Doublet
 CONVERGENCE THRESHOLDS:    1.00E-05 (Density)    1.00E-07 (Energy)
 MAX. NUMBER OF ITERATIONS:       60
 INTERPOLATION TYPE:            DIIS
 INTERPOLATION STEPS:              2 (START)      1 (STEP)
 LEVEL SHIFTS:                 -0.60 (CLOSED) -0.30 (OPEN)

 Number of closed-shell orbitals:    7 (   6   1 )

 Singly occupied orbitals:      8.1
 Orbital guess generated from atomic densities. Full valence occupancy:   10   2

 Molecular orbital dump at record        2100.2

 Initial alpha occupancy:   8   1
 Initial beta  occupancy:   6   1
 Wave function symmetry:    1

...

 ELECTRON ORBITALS
 =================

  Orbital  Occupation    Energy     Coefficients

...
    6.1     2.00000    -0.61166     1 2px   0.54607     2 2px   0.72110
    7.1     1.00000    -0.27385     1 2px   0.81431     2 2px  -0.59346
    8.1     1.00000    -0.15389     4 1s    0.63065     4 1s    0.39475

...
    1.2     2.00000    -0.58188     1 2py   0.54818     2 2py   0.67678
...

4) Eventually I found a workaround by invoking, after a first RHF calculation without occupation constrains, the MULTI program with a ROTATE card:

 {rhf; wf,15,1,1; orbital,2100.2; }

 {multi
 start,2100.2
 rotate 7.1, 8.1
 close,6,1
 occ,7,1
 wf,15,1,1
 orbital,2101.2}

This gives the desired result [[n2hrb+_rohf_multi_test.out]]

I will stick to the last solution but I was wondering what is the reason of the weird behaviours is the other cases, especially #3.

Best regards,

Alexandre

[Peterson, Kirk]

Dear Alexandre,

 

I'm not sure why this happens, but note you can also use the rotate command in the rhf program. So just follow your first simple rhf with another equally simple rhf with the rotate directive.

 

regards,

 

-Kirk

--
You received this message because you are subscribed to the Google Groups "molpro-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email to molpro-user...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/436ea328-9cc8-4de6-af23-c7a2d8183048n%40googlegroups.com.

[Alexandre Voute]

Dear Kirk,

Thank you for the notice, I missed the rotate directive in scf :)

I just tested your suggestion and it works indeed.

Best regards,

Alexandre