MRCI calculation of highly excited states of CO

[Gileadc]

Dear Molpro users and devs  

I am currently trying to calculate the potential energy curves of highly excited diatomic molecule states for use in Dissociative Recombination. I have seen a paper (by Steven L. Guberman,  "Potential Curves for the Dissociative Recombination of CO+", The Journal of Physical Chemistry A, 2013, 117, 9704−9711) which uses MCSCF and MRCI calculations using molpro to calculate the PEC of CO. Hence I am trying to replicate the paper's results to better understand how to implement the MRCI method but it seems to crash every time. For example, for calculating the PI triple states the paper says: "The orbitals were determined in complete active space self-consistent field calculations in which 1σ and 2σ were kept fully occupied and the active space consisted of the 3σ−6σ, 1π, and 2π orbitals" and "The CASSCF calculations are done in an average of 3 B1 and 3 B2 symmetry of C2v so that the πx and πy orbitals are equivalent. The orbitals are obtained from an eight-root average in the CASSCF. Eight roots are obtained in the MRCI in 3 B1 symmetry". It is also stated that the calculation uses the cc-pVQZ basis set.

therefore the MCSCF and MRCI inputs are:

***,co
MEMORY,100,M    
geometry={c;o,c,r}
basis=cc-pVQZ

dr=0.2
i = 1
do r=1.8,6,dr

{hf;wf,14,1,0;noenest;orbital,2100.2}

{MCSCF;start,2100.2;
Closed,2,0,0,0;
Occ,6,2,2,0;
wf,14,1,0;state,1
wf,14,2,2;state,8
wf,14,3,2;state,8
}
{MRCI;
wf,14,2,2;ref,3;state,8;
}

do j=1,8,1
e(8*(i-1)+j) = energy(j)
enddo

rr(i)=r
i = i + 1
enddo

{table,e
save,co_e_4.txt}

{table,rr
save,co_r_4.txt}
---

But it always seems to crash at the first radius - iteration in the MRCI program, giving me the following error:

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

CI VECTOR FOR STATE 6 DOES NOT OVERLAP SUFFICIENTLY WITH REFERENCE VECTORS

 It may help to increase NSTATI, e.g., use OPTION,NSTATI=9

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

I have tried to fix this by multiple ways but to no avail. Any suggestions? or would you recommend any other method for calculating highly excited states? 

Thank you in advance.

Regards,

Gilead

[paul.ad...@gmail.com]

I am attempting some similar calculations for N2, so I'm trying to debug your input deck and maybe we can learn together. I'm trying to run the input below and will let you know how it works. 

Why do you have 'wf,14,1,0;state,1' in the MCSCF block? 

Also, I'm starting the scan closer to the equilibrium distance, which seems to work better for me (then you just have to go in the opposite direction, too, to complete the scan). And I think you can just run HF once and then use previous MCSCF orbitals (the integrals will update) at later r values. I'm learning, too, for these excited state calculations, so the input below could very well have issues, too...we'll see.

====co_3Pi.inp====

***,co

MEMORY,100,M

geometry={c;o,c,r}

basis=cc-pVQZ

hf

dr=0.2

i = 1

do r=2.25,6,dr

{MCSCF

 Closed,2,0,0,0;

 Occ,6,2,2,0;

 wf,14,2,2;state,8

 wf,14,3,2;state,8

}

{MRCI;

 wf,14,2,2;ref,3;state,8;

}

do j=1,8,1

 e(8*(i-1)+j) = energy(j)

enddo

rr(i)=r

i = i + 1

enddo

{table,e

save,co_e_4.txt}

{table,rr

save,co_r_4.txt}

[paul.ad...@gmail.com]

Ah, I forgot to set r=2.25 before the hf line. It's running now.

[paul.ad...@gmail.com]

It looks like my results are very similar to what's in the paper for the 3Pi states. Thanks for posting this question...it was very instructive for me. Here is a plot (I only did r>2.25) and the input

Now, I have a question...can I get the natural orbital kinetic energy one-electron properties to print for the MRCI states?  When I try, it seems like Molpro is only giving me the kinetic energies for one symmetry.

====co_3Pi.inp===

***,co

MEMORY,100,M    

geometry={c;o,c,r}

basis=cc-pVQZ

dr=0.2

i = 1

r = 2.25

hf

do r=2.25,6,dr

{MCSCF

 Closed,2,0,0,0;

 Occ,6,2,2,0;

 wf,14,2,2;state,8

 wf,14,3,2;state,8

}

{MRCI;

 wf,14,2,2;ref,3;state,8;

}

do j=1,8,1

 e(8*(i-1)+j) = energy(j)

enddo

rr(i)=r

i = i + 1

enddo

{table,e

save,co_e_4.txt}

{table,rr

save,co_r_4.txt}

[Gileadc]

Thank you very much for your help!. I ran the calculation as well and got very similar results to the paper, as you said. 

As for your question, do you mean for this calculation or are you asking generally? (in this calculation we only calculate one symmetry in MRCI). If you are asking generally, can you please post an example code? 

[paul.ad...@gmail.com]

I was asking generally. I am interested in the NO KEs. I think I've figured it out by specifying `orbital,,natural,state=averaged` or `orbital,,natural,state=1`, etc in a property block with ekin.