CP-corrected interaction energy stays positive: what am I doing wrong?

[Joseph Van Forst]

Hi,

                                                 I’m trying to calculate the Rb–N₂ interaction energy at the CCSD(T) level in MOLPRO. For the basis I use N: aug-cc-pCVQZ and Rb: ECP28MDF with aug-cc-pVQZ-PP, plus a set of bond (ghost) functions placed between Rb and the N₂ center of mass. I calculate the energy of the full system at different Rb–N₂ separations, then perform the monomer calculations in the presence of the other basis functions. From these three calculations I build the counterpoise-corrected interaction curve:

E_int(r) = E(Rb+N₂) – E(N₂) – E(Rb)

The geometry is fixed (N–N = 1.1 Å, Rb along the axis, bond point at the midpoint), and I scan r = 3.4–7.8 Å.

The problem is that the curve is positive everywhere. It does show a shallow minimum around 6.5 Å, but the value at the minimum is still above zero. I would expect the well to be attractive (negative energy), as reported in other similar works. I also repeated the calculation without the midbond ghost functions, and the issue remains. 

My suspicions are:

1. Is eccsdt(i)=energt(1) the correct way to capture the total CCSD(T) energy in MOLPRO? Or should I instead use eccsdt(i)=energy right after {ccsd(t)} 

2. Is there any caveat in combining UCCSD(T) (for Rb+N₂ and Rb) with RHF-based CCSD(T) (for N₂) in counterpoise calculations?

3. Is there anything wrong with how I specify the wf cards? For example, for Rb I use wf,charge=0,symmetry=1,spin=1; uhf, and for N₂ I use wf,charge=0,symmetry=1,spin=0; hf. Could this be affecting the results?
   
   

Thanks for any advice. I'm adding all the input files to this mail. 
P.s: I tried the new approach with the "Interact" directive: but there I'm facing problem when adding the midbond function. That is the reason why I want to do it this way. 

[peterso...@gmail.com]

Hi,

 

Using energt(1) will work but just setting your variable to energy is a bit more standard.  Currently it doesn’t look like you’re correlating any electrons on Rb since you’re taking the default frozen core.  I would suggest using core,small.  I assume you also wanted to correlate the 1s of N since you’re using the aug-cc-pCVQZ basis set? 

 

You’re probably already aware, but there are correlation consistent basis sets available for Rb that are matched to the ECP28MDF pseudopotential. These are under the name aug-cc-pwCVQZ-PP, etc. in the Molpro library.

 

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, visit https://groups.google.com/d/msgid/molpro-user/fcebc02e-86b9-480a-92dd-7b9d1ff9f9c5n%40googlegroups.com.